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The dynamic analysis of large, complex structural systems is 
computationally intensive and therefore prohibits the use of optimization 
procedures, which are both iterative and complex with respect to variable 
search patterns. The solution to this problem is through the use of time and 
frequency synthesis techniques. They provide a means of rapidly recalculating 
a system's changed response due to structural modifications, as dictated by the 
optimization procedure. The efficiency is gained through the fact that the 
synthesis methods are independent of model size, in that only those model 
degrees of freedom where changes are made are required in the analysis. 
Furthermore, these methods are exact in their formulation, including the 
treatment of non-proportional damping. These structural synthesis 
techniques are developed in the context of optimal design of shock and 
vibration isolation systems. Their utility and value is demonstrated in the 
optimal design of an isolation system for a 109 dof non-proportionally 
damped structural system. In the course of the optimization, the synthesis 
techniques make possible 80 transient, frequency response, and static analyses 
in 2 hours and 39 minutes (desktop computer), while yielding an isolation 
design which satisfies all design constraints. 
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I. INTRODUCTION 



With the development of the finite element (FE) method, and the 
advancement of computer technology, it is now possible to design and 
analyze complex engineering systems. Through FE techniques, a detailed 
mathematical model of a complex structural dynamic system may be 
developed, and the static and dynamic responses determined. The 
'traditional' procedure for conducting a finite element analysis (FEA) of a 
structure is to first assemble the system matrices. These system matrices are 
the mass, stiffness, and damping matrices, which constitute the FE model of 
the complete structirre. The next step would be to determine the system 
responses using various solution techniques. This is the most 
computationally demanding phase of the FEA process. The results are then 
processed and interpreted. As a result of this analysis, the engineer may wish 
to change some aspect of the design and perform the analysis again. Even 
more useful would be the use of some type of optimization scheme to find 
the optimum design change while concurrently performing the FEA. 
However, if the 'traditional' method of FEA is used, every time a design 
variable is changed, the affected system matrices must be reassembled and the 
entire solution phase must be reaccomplished. Depending on the complexity 
of the system, and the number of different design parameters which may be 
changed, this route would be computationally impractical. As a result, the 
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designer may only be able to iterate through the design process a few times, 
and be left with a design which is less than optimum. 

Therefore, more efficient techniques for calculating a system's 
responses after modifications have been made, must be introduced. One such 
technique is the use of synthesis procedures to obtain the modified system 
responses. These synthesis procedures involve the use of the structure's 
baseline (pre-modification) responses, the modifications made to the 
structure's mass, stiffness, and damping matrices, and the equations which 
link the two to obtain the modified /synthesized responses. The advantage of 
this is that the system matrix assembly and solution phase of the FEA process 
must only be accomplished once. A smaller set of change matrices are 
assembled, and along with the presynthesized responses and computationally 
efficient synthesis equations, the synthesized responses are obtained. 

The two types of sythesis techniques to be featured in this thesis are 
frequency domain synthesis, and time domain synthesis. Frequency domain 
synthesis makes use of the baseline frequency response functions (FRFs) to 
calculate the new FRFs after the structure has been modified. Time domain 
synthesis uses the impulse response functions (IRFs) of the baseline structure, 
the coupling forces from the modification, and the convolution integral. The 
combination of these produces a nonstandard nonhomogeneous Vol terra 
integrodifferential equation (VIDE) of the second kind which is solved in 
order to calculate transient responses of the modified structure. The use of 
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these techniques greatly reduce the computer computation time and allow for 
the application of optinuzation techniques to complex engineering systems. 
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II. FINITE ELEMENT FORMULATIONS 



As mentioned previously, the initial step in the FEA process is to create 
a finite element model of the structure to be analyzed. In this study, different 
finite element types are used to ensure the validity and universal applicability 
of the synthesis techniques, and to assist in the synthesis formulations. For 
simplicity, linear lumped mass-spring systems were used. Its formulation is 
based on the use of Newton's laws to derive the equations of motion for each 
particle of mass [Ref. 1]. These systems allowed for the initial general 
formulation of the synthesis techniques and the building of a checking system 
to ensure its accuracy. In order to illustrate more realistic systems, systems 
formed from beam elements and plate elements were used. The beam 
element formulation is based on Euler-Bemoulli beam bending, and the use 
of Hermitian shape functions [Refs. 2 & 3]. The plate element formulation is 
based on shear deformation which includes both the transverse shear energy 
and the bending energy [Ref. 3]. One very important anomally about the plate 
shear deformation formulation is that, the unconstrained stiffness matrix (K) 
is rank deficient by the number of degrees of freedom per node, plus one. 
When solving the eigenvalue problem to obtain the system's natural 
frequencies (co^) and mode shapes, one additional "spurious" rigid body mode 
is obtained which disrupts all modal calculations. Therefore, to avoid this 
problem, all analyses done using the plate were done on a constrained to 
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ground system to eliminate all rigid body modes. The constraints were then 
subsequently synthesized out during the analysis in order to obtain the truly 
desired responses. 
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III. FREQUENCY DOMAIN STRUCTURAL SYNTHESIS 



The following background information on frequency domain 
structural synthesis is taken from references [4 & 5]. This portion of the thesis 
will outline the methodologies and formulations of the frequency synthesis 
technique and the classical techniques used to check its accuracy. 

Frequency domain techniques were demonstrated as early as 1939 by G. 
Kron in his book on the tensor analysis of electrical networks. Since then, 
numerous formulations and application techniques have been used and 
documented. The advantages in computational time, and the flexibility of 
this technique, have also been well documented [Ref. 6]. The solutions 
obtained through this method are exact with no approximations. 

As mentioned previously, frequency domain structural synthesis is a 
methodology whereby changes can be made to a given structure, and its new 
frequency response fvmction (FRF) can be formulated through the use of a 
single synthesis equation. This is quite different from the classical method of 
evaluating a structural change through the reconstruction of the basic 
elements which define the structure, and then using an impedance inversion 
or modal superposition technique on the new structure. The frequency 
synthesis technique may be divided into two major classes, coupling and 
modification. Coupling is the joining of a totally independent imcoupled 
structure to the given structure. Modification is the addition of redundant 
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load paths in the given structure. Each of the two classes of synthesis may be 
either direct or indirect. For the indirect case, there is the existence of an 
interconnection impedance element between substructures in a coupling 
operation, and degrees of freedom in a modification operation. In the direct 
case, the coupling and modification is done without the impedance element 
existing, therefore making the modification synonymous to applying a 
constraint. 

This thesis will focus on the indirect modifications to a given structure. 
It will allow for the inclusion of a spring stiffener or a viscous damper 
between two points, a lumped mass on the structure, and implementing a 
base excitation to the structure. A generalized computer program will be 
presented, which will allow for any of the above mentioned changes to be 
accomplished on some baseline structure. The program then returns the FRFs 
for all dofs where changes have occurred and any other dofs that may be of 
interest to the user. 

A. GENERALIZED FREQUENCY RESPONSE 

The initial step in obtaining the new synthesized FRFs for the modified 
structure, is to have the FRFs for the baseline structure. The following 
equations are for a system which is excited by a force on the structure or by a 
base excitation. These formulations are also used to verify the accuracy of the 
synthesis after a modification has been made. 
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Consider the following two degree of freedom system: 




F(t) 



Figure 3.1 Two Degree of Freedom Mass-Spring-Damper System 



where: 

mp = mass of a plate 

m^ = mass of computer equipment 

kp = isolator stiffness between plate and ground 

kj = isolator stiffness between plate and computer equipment 

Xp(t) = plate displacement as a function of time 

x^(t) = computer displacement as a function of time 

F(t) = excitation force 

Applying Newton's 2nd law, the equations of motion for this system may be 
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written as: 



mp 0 
0 



c J L 



Cp+^c -C, 

-c, c. 






k, -k 
k 



p 1-. 



C Jl.Xc. 



(3.1) 



Generalizing this formulation to an ndof (number of degrees of freedom) 
system where the external force may occur at any dof, and using a more 
compact matrix and vector notation, the 2nd order system of linear equations 
for an ndof structure can be written as: 



[M]{x} + [C]{x} + [K]{x} = {F} (3.2) 

Assuming a harmonic forcing function and consequently a harmonic steady- 
state response: 

{F(0} = {F}e^ {x(f)} = {x}e^ (3.3) 

where: 



{f} & |x] are the amplitudes of the harmonic forcing function and 



response. 

Taking the first and second derivatives of the response (x(t)}, and substituting 
into equation (3.1) and simplifying: 



[k + - q^m]{x} = {f} 


(3.4) 


[Z(n)]{x}={F} 


(3.5) 



[Z(Q)] is known as the impedance matrix. Multiplying both sides of the 
equation by [Z(Q)]‘^ and denoting this as frequency response function (FRF) 
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matrix [H(Q)], equation (3.5) becomes: 



{x} = [H(Q)]{f} 



(3.6) 



where: 



[H(Q)] = 



H21 H22 




(3.7) 




ndof X ndof 



and in elemental form 




5 

F; 



(3.8) 



The FRF matrix [H(n)] is both complex valued, and frequency dependent. The 
advantage of this formulation in this work is that it can be used to check the 



structure. These damper modifications represent non-proportional damping, 
which is not easily handled in modal coordinates. 

The modal coordinate formulation of the FRF matrix begins from the 
same general equations of motion presented in equation (3.2), and the 
assumptions of harmonic excitation and response in equation (3.3). It also 
includes the assumption of a harmonic modal response: 



frequency synthesis method when damper modifications are made to the 



{?«)}= {?k'“ 



(3.9) 



and the linear modal transformation: 



{x(r)} = [0]{^(r)} 



(3.10) 
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where [<I>] is the assembles matrix of mass normalized mode shapes. By 
applying equations (3.3), (3.9), & (3.10), to equation (3.2) and simplifying; 

[-Q"[M][<I>] + yS2[C][<l>] + [/iTlK]]® = {F} (3.11) 

Premultiplying equation (3.11) by [O]^: 



= (3.12) 





1 


+ jQ. 




+ 






_ 








_ * 



where it is recognized that; 

[<Df[M][0] = [I] 
[<]>f[C][0] = [2Ccol 
[Of[K][0] = [co'] 
[Of {F} = {5} 



(3.13) 

(3.14) 

(3.15) 

(3.16) 



^ is the modal damping factor, co is the natural frequency, and the subscript r 
represents each mode. Rewriting equation (3.12) and solving for the modal 
displacement yields: 



m= 



1 



0 ).^ - 

Using equations (3.10) & (3.16), equation (3.17) is transformed back to physical 



{ 5 } 



(3.17) 
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coordinates: 



{x}=[o] 



1 






[<I.f{F} 



(3.18) 



From equation (3.6): 



{H(Q)} = [<D] 



o>/ - 



[o>r 



(3.19) 



and any element of the FRF matrix may be represented as: 



nxnooo 

H„(n)= £ 



m 



(0, -ir +42C,co^Q. 



(3.20) 



Although this modal coordinate formulation does not handle non- 
proportional damping, it does provide for the ability to sum over a subset of 
the complete modes, rather than summing over all modes for the FRF. The 
number of modes necessary to correctly capture the response is dependent on 
the frequency range of interest. The lower the frequency range of interest, the 
less number of modes necessary. The obvious advantage of this procedure is 
that for large dof systems, this summation may be truncated, and 
computational time saved. However, the question of how many modes are 
enough must be considered and truncation criteria must be established. The 
modal formulation will be the one used for the s)mthesis techniques. 
Recognizing the modal truncation issue, this thesis uses all modes in its FRF 
calculations. 
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Now consider the following two degree of freedom system: 




Figiire 3.2 Base Excited Two Degree of Freedom Mass-Spring-Damper System 



where: 

kb = isolator stiffness between plate and base 
y(t) = base excitation displacement 

all other variables are the same as in the system illustrated in Figure 
3.1. 



Applying Newton's 2nd law, the equations of motion for this system may be 
written as: 



m 



P 0 

0 m. 



Cb+Cc -Cc 

-c, c. 



kb + k, -k, 

-k, k. 



PL- 



Cb 0 
0 0 



{y}- 



{y} 



kb 0 



0 0 _ 
(3.21) 



Generalizing this formulation to an ndof system where the base excitation 
may occur at any dof, and using a more compact matrix and vector notation, 
the 2nd order system of linear equations for an ndof structure can be written 
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as: 



[M]{ii} + [C]{x} + [K]{x} = [C,]{y} + [K,]{y} (3.22) 

Following the same procedure as with the FRF formulation for a force 
excitation on the structure, and assuming a harmonic base excitation 

({y(t)} = {y}c^), the following result is obtained: 

{X} = [k + + /QC,]{y} (3.23) 

Since the base excitation |y} is the same at all dofs, Y can be factored out and 
the result is: 



{X} = [K + J{1}Y (3.24) 

and 

{H} = [K + [K^ + ]{!} (3.25) 

In this instance the collection of FRFs is not a matrix but rather a vector 
where in elemental notation: 



H = 



Y 



(3.26) 



In some instances, the base excitation and desired response may not both be 
displacements. One different combination could be an interest in obtaining 
the output displacement response due to a given base acceleration. From the 
assumption of a harmonic base excitation and response: 

{x(t)} = {x}c^ 

{x(t)} = = -Q^{x}e^ 



{y(t)} = {Y}c^ 

{y(t)} = 
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(3.27a&b) 

(3.28a&b) 



(3.29a&b) 



The following table uses the above relationships between acceleration and 
displacement to illustrate possible combinations and how to handle each. 





Y 


Y 


X 


[H] 


[H] 




1 


X 




[H] 



Table 3.1 Response and Base Excitation Relationships 



The modal formulation for base excitation is not needed since the 
synthesis method for a base excitation or force excitation, uses the previous 
modal FRF formulation exclusively. 



B. FREQUENCY DOMAIN SYNTHESIS FORMULATION 

Now consider the following general ndof (number of degrees of 
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freedom) structure, to which structural modifications are to be made: 




Figiure 3.3 General NDOF Structural System 



The letters i, c, and b represent the physical coordinate system for the 
structure where: 

i = iset - the set of internal dofs where no changes occur, but there is an 
interest in knowing FRF information about these dof. 
c = cset - the set of dofs where changes have occurred, 
b s bset - the set of dofs where the structure is indirectly connected to a 
base excitation. 

Using equation (3.6) and the previously described coordinate system for 
rearrangement and partitioning, the structural system is described in the 
frequency domain at each frequency as: 



Xi' 



r[H„]J [Hj i [H,,]] 




[Hj ; [H„] ; [H,^] 




[[H«]l [H^] i [H jJ 





(3.30) 
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Before modification, the force vector refers to externally applied forces 
to the structure. However, following modification, there exists coupling 
forces between the cset and bset dofs. By definition, the iset dofs do not 
experience these coupling forces. Therefore, the partitioned forces may be 
rewritten as: 



fi=^“‘ (3.31a) 

f^ = f + f; = f“‘ - AZ,x; (3.31b) 

k = r + k = r - AZ,(x; - y) (3.31c) 



where: 

P - the coupling forces due to structure modification. 

[AZ] - the impedance matrix which is equal to [AK + jQAC + Q^AM]. 
x’ - generalized (synthesized) responses after modification. 
Substituting equations (3.31) into equation (3.30), and recognizing that the 
presynthesized responses {xj x^ Xj,}^ are those due to the external forces only, 
the following result is obtained: 



bj 



X; 



> =<x. 



[H.] [H»] 
[H„] [H^] 
[[Hb,] [Hbbll 



[AZ,] 0 ■ 

0 [AZ,] 



[HJ [Hi,] 
[H„] [H„] 
[[Hb.1 [H„]> 



[AZ,] 0 ■ 

0 [AZ,] 



fO] 

yj 



I: 



(3.32a) 



Denoting two additional coordinate sets 'e' and 'z' where: 



e = eset = iu cu b 



z = zset = cu b 
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and [A^ is the total impedance matrix. Equation (3.32a) can be rewritten as: 



KF = +[H.jAZ»]{y} (3.32b) 

From equation (3.32b) it can be seen that there are two imknown synthesized 
responses in this one equation. Another independent equation can be 
generated from equation (3.32a) without introducing any additional 
unknowns. By observing the bottom two rows of equation (3.32a), the 
following relationship is derived: 

{x.}' = {xJ-[H„][AH]{xJ- +[H.jAZ,]{y} (3.33) 

Solving for {x^}* and substituting that expression into equation (3.32b) yields 
the general expression for the responses of the synthesized structure. 

{x.F = {x.} - [H„][AZ]([I + H„AZ]-*{x J +[I + H„AZ]-'[H*][AZ J{y}) + [H.J[AZ,]{y} 

(3.34) 

From the general equation for frequency response, it is recognized that: 

{x.}=[H„]{fr} {x,} = [H„]{r} (3.35a&b) 

Equation (34) is therefore rewritten as: 

M = [H„]{r} - [H«][AZ]([I + H^AZ]-’[H„]{r} + [I + H^AH]-’[H^][AZ,]{y}) + [H.,][AZ,]{y} 

(3.36) 

As can be seen from equation (3.36), the synthesized response is now in terms 
of the known FRFs of the presynthesized structure, known impedance due to 
modifications, and known force and base excitations. 
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Now that the complete general frequency response synthesis equation 
has been derived, consider the case where there is no base excitation applied 
to the structure. By letting {y} = 0, the general s)mthesis equation becomes: 

K}' =([H„]-[H„lAZ][l + H„A2]-'[H„]){fr} (3.37) 

Since no base excitation exists, the impedance between the structure and base 
(AZb), is reduced to a V dof to groimd change. Therefore the bset dofs are 
really additions to the cset: z = cub = cuc = c and AZj, = AZ^. Also 

recognizing that {xg}* =[Hgg]*|fg““j, the synthesized FRF matrix at each 
frequency can be written as: 

[H.]' =[H„]-[H^][AZjI+H,,AZj-'[H„] (3.38) 

Now consider the case where there are no external forces applied to the 
structure. By letting {fe**‘} = 0, the general synthesis equation becomes: 

{x.}' = ([H,jAZj-[H„][AZII + H^4Z]-'[H,^][AZj){y} (3.39) 

Factoring out the base excitation magnitude y and applying equation (3.6) yet 
again: 

{H,}'=([H.J[AZ,]-[H„][AHII + H„A2]-'[H^][AZ,]){1} (3.40) 

In elemental notation, the FRF vector represents the following at each 
frequency: 

{H.r={x.}-i (3.41) 
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Equations (3.38) and (3.40) show how the new synthesized FRFs are 
obtained from the presynthesized FRFs and the impedance caused by 
structural modifications. These synthesis equations are exact, and there are no 
limitations to the modification values. The change may even be negative to 
represent the removal of mass or removal of an interconnection element 
from the structure. 

C FREQUENCY DOMAIN SYNTHESIS COMPUTER CODE 

The computer language used for the frequency domain synthesis and 
all other computer coding in this thesis is MATLAB V.4.2c.l. The goal in the 
formulation of the frequency domain synthesis computer programs is to 
perform the synthesis for an indirect modification operation to a given 
structure. The program allows for the inclusion of a spring stiffener or a 
viscous damper between two dofs, a lumped mass addition on the structure, 
and implementing a base excitation to the structure. The program then 
returns the FRFs for all dofs where changes have occurred and any other dofs 
that may be of interest to the user. Appendix A shows the computer codes 
used to perform the frequency domain synthesis, and perform a comparative 
analysis of the synthesis versus classical methods. The term classical is 
synonymous to the traditional method of reformulating the elemental 
matrices following a modification. 
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The programs presented in Appendix A can be segregated into three 
categories. The first is the calculation of the synthesized FRFs due to an 
external force excitation on the structure, and their comparison to the 
classically calculated FRFs. The second is the calculation of the synthesized 
FRFs due to a base excitation on the structure, and their comparison to the 
classically calculated FRFs. Both of these programs use the impedance 
inversion method to calculate the presynthesized FRFs, and the classical FRFs 
following modification. The third category is the calculation of the 
synthesized FRFs where the different aspects of the synthesis technique have 
been modularized into MATLAB functions. This modularization assists the 
process in running more efficiently, and allows for the universal use of the 
frequency synthesis technique. This is crucial since the ultimate goal of the 
technique is to be used efficiently in another process (i.e. optimization). Since 
efficiency is of concern, the modal method of calculating the presynthesized 
FRFs is used. 

The following is a diagram of the flowpath of the modularized 
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frequency synthesis programs: 



LoadFRFs 

changejn 



Load Structure 


-> 


Calculate FRFs 




__aE 

Modifications 




Form Change 
Matrices 


— > 


Perform Synthesis 




FRF Selection 


changejn 




modaLm 

frftnodaljn 




changejn 




delta jn 




frfsynthjn 




frf_siftjn 



Figure 3.4 Flowpath of the Modularized Frequency Synthesis Programs 



A more detailed description of each function is contained within Appendix 
A. 

D. ILLUSTRATIVE EXAMPLES 

The purpose of this section is to illustrate the use of the frequency 
domain synthesis technique and validate its accuracy. The computational 
efficiency of this technique has already been well documented [Ref. 6], 
therefore a time comparison is not done. The three structural systems 
previously described, mass-spring-damper, beam and plate, are used to 
accomplish this goal. 
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1. Mass-Spring-Damper System 

Consider the following mass-spring-damper system: 




k = 1000 Ib/in 
= 250 IbAm 




^ k 






m 




m 



m = 0.259 Ib-sVin 
C30 = 2 Ib-s/in 



x(t) 



Figure 3.5 Mass-Spring-Damp>er System Experiencing Force Excitation 

The solid elements represent the original structure and the dashed elements 
represent modifications. The subscript 'i' in the excitation force represents the 
dof where the excitation is applied. Although not shown in the figure, the 
original structure is proportionally damped using [C] = a[K]. The 
modifications are arbitrarily selected and the addition of the damper element 
represents nonproportional damping. The following is the resultant 
synthesized and classically determined FRF which represents the 

response magnitude at dof 3 due to a force at dof 2. 
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Figture 3.6 Force Excitation Mass-Spring-Damper System FRF 



From the graphs produced, it can been seen that the plots are identical, 
proving that the synthesis technique is exact, and that the computer coding 
for the arbitrary changes is correct. These results were produced using the 
fsynstr.m program. 
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2. Beam System 

Now consider the following beam system: 




7 

/ / / 



= 1000 Ib/in 
Cj 3 o = 20 Ib-s/in 
m„ = .5 Ib-sVin 



Figure 3.7 Beam System Experiencing Base Excitation 



where the beam parameters are : 

Length: L = 5 ft Width: w = 3 in Height: h = 4 in 

Young^s Modulus: E = 10e6 psi 
density: p = 2.53e-4 lbf-sec'^2/in^4 

The free-free beam which is discretized into 10 beam elements, 11 nodes, and 
22 dofs, represents the original structure. The modifications are represented 
by the dashed lines. Just as in the previous example, proportional damping is 
used to form the original [C] matrix and all changes to the beam are arbitrary. 
The following is the resultant synthesized and classically determined FRF 
Hj 3 (Q) which represents the response magnitude at dof 13 due to the base 
excitations. 
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Frequency Response Function (HI 3) 




Figure 3.8 Base Excitation Beam System FRF 



Just as before, it can be seen that the two plots are identical. These results were 
produced using the fs)mbase.m program. 



3. Mass-Plate System 

Up until now the structures use to demonstrate the frequency synthesis 
technique were relatively simple with relatively small number of dofs. The 
true test of the technique and the computer code is in their applicability to 
relatively large structures. Therefore, consider the following computer-plate 
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system: 





Cb 






y(t) 



J 



Figure 3.9 Plate-Mass System Experiencing Base Excitation 



where the computer and plate parameters are : 

Plate: Width: w = 5 ft Depth: d = 5 ft Thickness: t = 1 in 

Young's Modulus: E = 30e6 psi Poisson's Ratio: v = 0.3 

Density: p = 7.35e-4 lbf-sec'^2/in^4 
Computer: Weight: W = 100 lbs 

In this case, the original structure is represented by the free-free plate 
with the attached computer located off center and above the plate. The plate is 
discretized into 25 plate elements, 36 nodes, and 108 dofs. The computer is 
modeled as a single lumped mass in the translational direction only, and 
increases the total number of nodes and dofs to 37 and 109 respectively. The 
modifications are represented by the dashed lines and are comprised of 
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changes to the isolators between the base and the plate. Just as in the previous 
example, proportional damping is used to form the original [C] matrix, and 
any changes to the system are arbitrary. However, since the manipulation of 
this system in an optimization routine is the ultimate goal of this research, 
the isolator elements which will act as optimization variables are chosen for 
modification. Therefore, modifications were also performed on the isolators 
between the plate and the computer. The modularized program which uses 
the modal method of calculating the original FRF will be used in this 
instance. However, the classical method will also be used as a check system. 
The following is the resultant synthesized and classically determined FRF 
H^o^CQ) which represents the response magnitude at dof 109 (the computer) 
due to the base excitations. 
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Figure 3.10 Free-Free Plate-Mass System FRF (Modal Synthesis) 

From the plots it is easily observed that the synthesis method and the 
classical method do not yield the same results. After successfully performing 
countless other frequency synthesis analysis on all three structures, and 
successfully re-rimning this same analysis where the original FRF is obtained 
using impedance inversion methods vice modal methods (Fig. 3.11 shows). 
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Figure 3.11. Free-Free Plate-Mass System FRF (Impedance Synthesis) 

it is concluded that the problem lies in the use of the modal method of 
calculating the original FRF. After observing the mode shapes and natural 
frequencies of the original free-free plate-mass system, the existence of an 
additional spurious rigid body models discovered. The existence of this rigid 
body mode vice a flexible mode makes the FRF calculations incorrect. After 
unsuccessful attempts at trying to replace the rigid body modes naturally 
generated by this finite element formulation, with rigid body modes created 
using the Graham-Schmidt method [Ref. 2], the decision to constrain the 
plate-mass and then remove the constraints using synthesis is reached. What 
this accomplishes is to eliminate the rigid body modes entirely and allows for 
the proper calculation of the original FRF. However, there is a computational 
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price to pay with the addition of four more load paths to contend with in the 
calculations. This is however, a very minute price to pay when compared to 
using the impedance inversion method. 

The following is the resultant modally-synthesized and classically 
determined FRF when the original plate-mass structure is exactly as 

before, but with 1000 Ib/in spring-to-groimd elements located at all four 
comers. In the process of the analysis, these elements are synthesized out of 
the structure. Other modification values are still the same as before. 




Figure 3.12 Corners-to-Ground Plate-Mass System FRF (Modal Synthesis) 



It can be seen from these plots that eliminating the rigid body modes in the 
original structure, and then synthesizing the constraint elements out is an 
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accurate means of avoiding the problem with the rigid body modes. Another 
verification of the ability to synthesize out elements is that fact that these 
plots also match the plots in Figure 3.11 (free-free plate-mass system); as well 
they should. This will therefore be the method used in order to handle the 
dynamic analysis of the plate-mass system. 
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IV. STATIC DISPLACEMENT SYNTHESIS 



When using optimization procedures for the design of a system, there 
are usually various aspects of the design problem which are at odds with one 
another. For example, in the design of a support beam, the aspect to be 
optimized may be the weight of the beam. One method to reduce weight 
would be to reduce the beam's cross sectional area. However, since this is a 
support beam and will therefore be carrying some sort of load, a reduction in 
cross section results in an increase in stress. Since there are material 
limitations on the amovmt of stress the beam can experience, the reduction in 
cross sectional area is constrained by the stress limit. 

For a shock and vibration isolation systan, one aspect of the design 
which should be considered is the static sag or displacement of the system. 

The static sag could be used in the optimization problem as a constraint 
which limits the amount the isolator stiffness may be modified. Since the 
solution for the static response of a structure involves the inverse of the 
stiffness matrix ([K]"’), for large structures, it would not be computational 
efficient to perform this operation every time the optimization variables are 
changed. Therefore, some sort of synthesis technique must be employed to 
alleviate the need to calculate [K]'’ directly. This portion of the thesis explores 
the methodologies and formulations for the use of frequency domain 
synthesis in the calculation of static displacement. From henceforth, this 
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technique will be referred to as static displacement synthesis. Classical 
techniques such as Guyan model reduction are used to check the accuracy of 
these formulations. 

A. GENERALIZED STATIC DISPLACEMENT 

From Newton's 2nd law, the static equations for an ndof system is: 

{F} = [K]{x} (4.1) 

By defining the force vector as the weight of the structure, and assuming that 
there are no external forces on the structure, the displacement vector {x} 
represents the static displacements of the structure due to its own weight. 
Therefore, the static displacements of the structure due to its own weight 
equals: 

(xU = [K]-’(F) = [K]-'[Ml(g) (4.2) 

where {g} represents a vector where the gravitational constant appears at all 
dofs that are affected by gravity (i.e. vertical translational dofs). 

Equation (4.2) returns the static displacements at all dofs of the 
structure. The size of the stiffness matrix of which the inverse is taken is ndof 
by ndof. Since the repetitive calculation of [K]*’ is unsatisfactory due to the 
possibility of its large size, and a relatively small subset of the static 
displacements of the structure are desired, the use of model reduction and 
static condensation schemes are investigated. 
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B. GUYAN REDUCTION / STATIC CONDENSATION 

The following concepts and equations for model reduction are taken 
from references [2 & 7]. The Guy an reduction defines a transformation matrix 
[T] which can operate on a subset of the structure's displacements and return 
the full ndof set of displacements for the structure. 



a = aset - those dofs arbitrarily selected as the set of active dofs. 
o = oset - the remainder of the dofs omitted from the aset. 

Substituting equation (4.3) into equation (4.1), and premultiplying both sides 
of the equation by [T]^ yields: 




(4.3) 



where: 




(4.4) 



{F} = [K]{x J 



(4.5) 



where: 



[K] = [Tf[K][T] 
{F} = [Tf{F} 



(4.6a) 

(4.7a) 



Therefore, the static displacements of the desired dofs equal: 



{xJ=[KrMF} 



(4.8) 
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It can be seen from equation (4.8) that the main cost of obtaining the desired 
static responses of the system is now the inverse of the reduced stiffness 

matrix [K] which has size of aset by aset vice the inverse of the full stiffness 
matrix which has size ndof by ndof. There is a definite improvement in 
computational efficiency, but the issue of forming the transformation matrix 
[T], which contains the inverse of the sub-matrix [K„J still exists. For a large 
structure, the size of the aset will usually be substantially less than the size of 
the oset, thereby making [K„„] ’ a time demanding operation. In an 
optimization procedure, [K^o] ’ can be performed prior to the iterations begin, 
and passed into the iteration portion of the optimization program. Even 
though [KjJ ’ must only be performed once, the formation of [T] is still 
imdesirable due to the possibility that FRF data is readily accessible and the [K] 
matrix is unavailable. Therefore, a method of obtaining the reduced stiffness 
matrix and reduced force vector from FRF data is desired. 

C STATIC DISPLACEMENT SYNTHESIS FORMULATION 

During the derivation of the frequency domain synthesis, the eset 
coordinate set was defined as iucub. This means that the eset includes the set 
of all dofs that the user is interested in, dofs where modifications occur, and 
dofs where a base excitation is attached. In other words, when looking at the 
structure as a whole, eset is synonymous to the set of active dofs or aset. 
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1. Formulation of Reduced Stiffness Matrix [K*] 

Using the aset/ oset coordinate system previously described but with 
the eset equaling aset, the expanded version of equation (4.6a) is written as: 

[K] = [k^-K^K^-X] (4.6b) 

Now using the same eset/ oset coordinate system and the fact that [H] = [Z]'^ 
the following relationship is written: 





[[Z-] j [z»]l 


M 


0 


[[h»] : [H..U 


[[z»] 1 [z~lj 


0 


W. 



Carrying out the matrix multiplication of equation (4.9), the first row yields 
the following two relationships: 

HeeZee+H^Zoe = I H^Zeo + H^Zoo = 0 (4.10a&b) 

Solving for in equation (4.10b) and substituting it into equation (4.10a) 
yields: 

H.Z.-H„Z„ZJZ„=I (4.11) 

Solving for yields: 

[H.]=[Z.-Z„ZjZ«]'' (4.12) 

Recognizing that Z = K +^C - Q^M, and applying the static condition that 
Q=0: 

[H„(0)] = [K„-K„K„-'K^f (4.13) 
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From equations (4.6b) & (4.13), it is determined that: 

[K] = [H^(0)]"’ (4.14a) 

Therefore, since the reduced stiffness matrix can be determined from an FRF 
matrix, when a structure is modified, the new reduced stiffness matrix can be 
determined from a synthesized FRF matrix. 

[K*] = [h:,( 0)]'’ (4.14b) 



2. Formulation of Reduced Force Vector {F*} 

Recognizing that (F) = [M]{l}g the expanded version of equation (4.7a) 
in eset/ oset partitioning is: 

(F) = [M„ - i M„ - K„K^-'M^]{l}g (4.7b) 

Another very important concept of Guyan reduction is the reduced mass 
matrix. 



[M] = [Tf[M][T] 

= [M^ - - K^K^-’M^] 






By comparing equations (4.7b) and (4.15b), it can be seen that: 

{F} = [M]{lJg 

if and only if 






W = [T]{lJ= ; 



(4.15a) 

(4.15b) 



( 4 . 16 ) 



( 4 . 17 ) 
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where {1^} is an eset length vector of all ones, and {l^} is an oset length vector 
of ones and zeros, with ones at the translation dofs and zeros at the rotations. 
In effect, what equation (4.17) is stating is that for every row of the 
transformation matrix, the sum of its columns is exactly one or zero. In 
general, this occurrence is not true. It is very dependent on which dofs are 
chosen as the eset. Therefore, it must be determined which dofs are allowable 
choices for the eset in order for this formulation to work. 

Consider the general stiffness matrix of a restrained structure: 

(4.18a) 



where represents the stiffness matrix of the unconstrained structure and 
represents the grounded stiffness matrix. Equation (4.18a ) can be rewritten 
in eset/oset partitioning format as: 



[K"] = 





K] 


1 

1 

1 


K] 










Kl 


1 

1 

1 








0 ; [K= ]J 



(4.18b) 



The zeros in the grounded matrix are due to the fact that grounds produce no 
off diagonal information. Now a displacement vector is chosen which 
satisfies the static equations of equilibrium where only forces appear in the 
active or eset dofs, and which produces no strain energy internal to the 
restrained structure. 
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From equation (4.19) it can be seen that: 




(4.20a) 



The vector which satisfies equation (4.19) is already chosen by default. The 
vector equals the previously defined ||*|- Therefore, {'FJ = {1^} and 



W = iV- 



Now consider the following portion of equation (4.19); 




(4.19) 



Repartitioning this equation into translational and rotational dofs: 



'Kl 


0 


0 


0 


[{'I'.ltl 






0 


Ki 


0 


0 


, {'j'.Jr 






0 


0 


Ki 


0 


' {'•'oIt 




0 


0 


0 


0 




[KU 




0 



where: 







IMr] 


KIk 




Wr 






Wx 


({■J'Jr. 




[wJ 



Case I - Taking the third row of equation (4.21) yields: 



(4.21) 



(4.22) 



Since [K= ] is a diagonal matrix, the only way equation (4.22) is true is if 
=[0]* Therefore, this proves that there cannot be any translational 
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grounds in the oset. Consequently, this means that all translational grounds 
must be in the eset. 



Case n - Taking the fourth row of equation (4.21) yields: 

[K»kR}. = [K=L{0} = {0} (4.24) 

Equation (4.24) is always true regardless of Therefore, this shows that 

there are no restrictions on rotational groimds being in the oset. Thus far it 
has been determined what must be in the eset, all translational grounds . It 
must now be determined what else is allowed in the eset. 



Consider the zero internal strain energy requirement: 



K!i 



kY 



14^. 



[0| 

lo| 



(4.20b) 



Taking the first row of equation (4.20b) following matrix multiplication 
yields: 

[K"]{'P.} = -[K“]{>rJ <4.25a) 

Repartitioning this equation into translational and rotational dofs: 



K"' 



I 
I 



kY 

“JIB 






RRj 



[{■J'.ItL 

[{'p.IkJ' 





1 

1 


K1 




xrj 




[Kli] 


\ 

1 

RT ! 


K] 



IR, 

RR. 



[{'•'.ItI 

[{'•'oIrJ 



(4.25b) 



Case in - if eset includes every translational dof and nfi rotational dofs, 
equation (4.25b) becomes: 



F-L 0 
0 0 



[Mr ■■ 



0 [K~L' 



0 



0 



.{0]J 



(4.26) 
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After matrix multiplication, equation (4.26) requires that: 



= {0} M-27) 

which means that for every row, the sum of the translational dof columns for 
an unrestrained structure equals zero. This condition is always true. 

Therefore, this shows that all translational dofs are allowed to be included in 
the eset. 

Case IV - if eset includes every translational dof and one rotational dof, 
equation (4.25b) becomes: 



kl 



K^1 

Kl] 

“ JRRJ 



Ri-ItL 

[WrJ 



K“ 



TR 

RR, 



'WtI 



( 4 . 28 ) 



After matrix multiplication and the fact that [K^]^{le}j ={0}, equation (4.28) 
requires that: 



Since there is only one rotation in the eset, the dimensions of [K^]^ are the 
# of eset translations by 1. Consequently, in order for equation (4.29) to be 
true, must always equal zero. By the definition that this case includes a 

rotational dof value in the eset, ^ [0], Therefore, one rotational dof 

may not be included in the eset. 

Case V - if eset includes every translational dof and greater than one 
rotational dof, the same requirement (equation 4.29) as stated in Case IV 

exists. However, for this case, does not have to identically equal zero. 
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but the sum of the columns for every row of [K^]^ must equal zero. In 

general, this sum is not zero. Therefore, it is concluded that no rotational dofs 
can be included in the eset. 

Case VI - if eset includes a subset of all translational dofs, equation 
(4.25b) becomes: 



'Kin 


o' 


r 


0 


0_ 


Ui.}kJ 1 



Kin [K~L 



0 



[{0} I 



(4.30) 



After matrix multiplication and rearrangement, equation (4.30) requires that: 



[K-L{I-}t+KLWt = {0} («1) 

Equation (4.31) is the same as if summing all translational columns of an 
unrestrained structure for every dof in the eset. From Case in (equation 4.27), 
it is seen that this sum is zero and equation (4.31) is satisfied. Therefore, this 
shows that any subset of all the translational dofs may be included in the eset. 

As a result of Cases I- VI, the following requirements on the choices of 
eset dofs for static displacement synthesis are summarized: 

(1) eset must include all translational groimded dofs 

(2) structure may be grounded an)rwhere 

(3) any subset of, or all translational dofs may be included in the eset 

(4) nii rotational dofs are allowed in the eset 

For the work to be presented in this thesis, the requirements necessary for the 
reduced force vector to be determined from the reduced mass matrix are met. 
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Therefore, the task is now to derive a methodology for extracting [M] from 
FRF data. 



3. Formulation of Reduced Mass Matrix [M*J 
Since [H^(0)] = [K]"\ it follows that: 



[H^(n)] = [K-Q'Mr 

Taking two derivatives of equation (4.32) with respect to Q yields: 



(4.32) 



d"H^(n) 

dQ" 






dH. 






= [K - [2QM] + [iC - a^M]-' [2M] + [20M] [K - 



dQ 



da 



(4.33) 



Applying the static condition of Q=0 to equation (4.33), and solving for [M]: 

[M] = i[a^^^51[K] (4.34a) 

2 dQ^ 

or rewritten in terms of the synthesis process: 

[M'l = |[H:.(0)f’^^5>[H:.(0)]-' (4.34b) 

From equations (4.34a&b), the unknown second derivative of the eset FRF at 

zero frequency is present. The calculation of this value will be handled using 

the following forward finite differencing scheme with error order O(Aa^): 

d^H:,(o) -H:.(n,)+4H:.(n,)-5H:.(a,)+2H:.(o) , 

dn" {anf 

The accuracy of this calculation is very sensitive to the size of AQ. The proper 
choice of AQ. is dependent on the relative magnitudes of the mass and 
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stiffness matrices. It was found for this work that AQ = 0.1 provided a very 

d^H* (0) 

accurate estimation of ^ — . However, this value is not universal and the 

proper selection of AQ should be based on individual structures. 

From the previously discussed requirements on the conversion of [M] 

to {F}, it can be seen that static displacement s)mthesis is not arbitrary in the 
selection of which dofs may be included in the eset. However, this restriction 
does not apply to which dofs can be chosen to be modified, but only to 
whether static displacement information can be obtained about these dofs. For 
example, modifications to rotational dofe may be made to a structure and the 
synthesized FRF obtained. As a result of this, rotational dofs exist in the eset, 
which is not allowed. The solution to this dilemma is to re-synthesize the 
structure with no changes and only translational dofs in the iset. What this 
does is remove the rotational dofs from the eset and condense the FRFs prior 
to using them to obtain the reduced stiffness matrix and force vector. 
Therefore, the static displacement method is really not restrictive in its 
application, but rather in what information it can provide. 

D. STATIC DISPLACEMENT SYNTHESIS COMPUTER CODE 

The goal in the formulation of the static displacement synthesis 
computer programs is to first perform the frequency domain synthesis for an 
indirect modification operation to a given structure. The program allows for 
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the same type of modifications to the structure that were presented in the 
frequency domain synthesis chapter of this thesis. The program then uses the 
synthesized FRFs and returns the static displacements for all dofs designated 
in the eset. 

The reason why the frequency domain synthesis is performed here 
seperately vice using the FRFs obtatined from doing a frequency domain 
synthesis problem is two fold. First of all, there is no damping effect for a 
static problem. Therefore, the synthesized FRFs must be generated with zero 
damping. Secondly, if we want to determine the static displacement for a base 
excitation problem, the base excitation form of frequency domain synthesis is 
not what is desired. The correct form of the frequency domain synthesis is the 
force excitation synthesis equation, with the base excitation interconnection 
elements considered as springs to ground. Furthermore, an entire spectrum 
of FRF information is not necessary for static displacement calculations. 
Primarily only the static condition of Q=0 is needed. However, in this 
implementation, the synthesized FRFs at the first four frequencies are needed 



for the finite differencing calculation of 



d^H:,(0) 



Appendix B shows the computer codes used to perform the static 
displacement synthesis, and perform a comparative analysis of the synthesis 
versus classical Guyan reduction methods. The programs presented in 
Appendix B can be seperated into two categories. The first category is a 
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program which uses the classical Guyan reduction techniques to calculate 
static displacement. The second category is the calculation of the synthesized 
static displacements where the different aspects of the synthesis technique 
have been modularized into MATLAB fimctions. The following is a diagram 
of the flowpath of the modularized static displacement synthesis programs; 




Figure 4.1 Flowpath of the Modularized Static Displacement Synthesis 

Programs 



It can be seen that the flowpaths for static displacement and frequency domain 
synthesis are very similar. This is done intentionally so that all synthesis 
techniques can be driven by the same modifications to a structure. A more 
detailed description of each function is contained within Appendix B. 



E ILLUSTRATIVE EXAMPLES 

The purpose of this section is to illustrate the use of the static 
displacement synthesis technique and validate its accuracy. The mass-plate 
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structural system used in the frequency domain synthesis example, along 
with the changes that were made to the structure, is also used here to 
illustrate static displacement synthesis. The following tables are comparisons 
of the reduced matrices, vectors, and static displacements calculated using 
frequency synthesis and Guyan reduction. 
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l . Oe+04 * 



[k*]=[h:,(o)] 



1.1518 


- 0.0000 


- 0.0000 


1.0000 


- 0.0000 


- 0.0000 


0.2276 


- 0.0000 


- 0.3794 


- 0.0000 


0.0000 


- 0.0000 



- 0.0000 


0.2276 


- 0.0000 


- 0.0000 


1.0000 


- 0.0000 


- 0.0000 


1.3415 


- 0.0000 


- 0.5691 


0.0000 


0.0000 



- 0.3794 


0.0000 


- 0.0000 


0.0000 


- 0.0000 


0.0000 


- 0.5691 


0.0000 


1.4985 


- 0.5500 


- 0.5500 


0.5500 



Table 4.1a Static Displacement Synthesis Reduced Stiffness Matrix 



[K] = [Tf[K][T] 



1 . 0e - K)4 * 




1.1518 


0.0000 


0.0000 


1.0000 


- 0.0000 


- 0.0000 


0.2276 


0.0000 


- 0.3794 


- 0.0000 


0 


0 



- 0.0000 


0.2276 


- 0.0000 


0.0000 


1.0000 


- 0.0000 


0.0000 


1.3415 


- 0.0000 


- 0.5691 


0 


0 



- 0.3794 


0 


- 0.0000 


0 


0.0000 


0 


- 0.5691 


0 


1.4985 


- 0.5500 


- 0.5500 


0.5500 



Table 4.1b Guyan Reduction Reduced Stiffness Matrix 
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d"H* (0) -H* (a)+ 4H* (fl») - 5H* (il) + 2H* (0) 



dQ" 



l.Oe-07 * 






0.0671 


0.0294 


0.0294 


0.0294 


0.0588 


0.0147 


0.0294 


0.0147 


0.0588 


0.0271 


0.0294 


0.0294 


0.1094 


0.0593 


0.0593 


0.1470 


0.0593 


0.0593 






0.0271 


0.1094 


0.1470 


0.0294 


0.0593 


0.0593 


0.0294 


0.0593 


0.0593 


0.0774 


0.1414 


0.1978 


0.1414 


0.3567 


0.5049 


0.1978 


0.5049 


0.8242 



Table 4.2a Static Displacement Synthesis 2nd Derivative Matrix 



d"H^(0) 

d£^ 



= [K][2M][K] 



l.Oe-07 * 




0.0671 


0.0294 


0.0294 


0.0588 


0.0294 


0.0147 


0.0271 


0.0294 


0.1094 


0.0593 


0.1470 


0.0593 



0.0294 


0.0271 


0.0147 


0.0294 


0.0588 


0.0294 


0.0294 


0.0774 


0.0593 


0.1414 


0.0593 


0.1978 



0.1094 


0.1470 


0.0593 


0.0593 


0.0593 


0.0593 


0.1414 


0.1978 


0.3568 


0.5049 


0.5049 


0.8242 



Table 4.2b Guyan Reduction 2nd Derivative Matrix 
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[M-] = |K(o)]" [h:.(o)]'’ 



0.1928 


0.0903 


0.0903 


- 0.0492 


0.0422 


- 0.0000 


0.0903 


0.2940 


0.0735 


0.0620 


0.1417 


- 0.0000 


0.0903 


0.0735 


0.2940 


0.0620 


0.1417 


- 0.0000 


- 0.0492 


0.0620 


0.0620 


0.1537 


- 0.0095 


- 0.0000 


0.0422 


0.1417 


0.1417 


- 0.0095 


0.4215 


- 0.0000 


- 0.0000 


- 0.0000 


- 0.0000 


- 0.0000 


- 0.0000 


0.2588 



Table 4 . 3 a Static Displacement Synthesis Reduced Mass Matrix 



[M] = [Tf[M][T] 



0.1928 


0.0903 


0.0903 


- 0.0492 


0.0422 


0 


0.0903 


0.2940 


0.0735 


0.0620 


0.1417 


0 


0.0903 


0.0735 


0.2940 


0.0620 


0.1417 


0 


- 0.0492 


0.0620 


0.0620 


0.1537 


- 0.0095 


0 


0.0422 


0.1417 


0.1417 


- 0.0095 


0.4215 


0 


0 


0 


-0 


0 


0 


0.2588 



Table 4 . 3 b Guyan Reduction Reduced Mass Matrix 
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SjTithesis: {F*} = [M*]{lg}g 



Guyan: 



{F} = [Tf{F} 



Synthesis: 



-141.6(XK) 


-141.6012 


-255.6080 


-255.6102 


-255.6080 


-255.6102 


-84.5960 


-84.5966 


-285.0197 


-285.0226 


-99.9984 


-100.0000 



Table 4.4 Reduced Force Vectors 



L=[K-r{F-) 


Guyan: {xJ^^,=[K] 


-0.0296 


-0.0296 


-0.0256 


-0.0256 


-0.0256 


-0.0256 


-0.0316 


-0.0316 


-0.0714 


-0.0714 


-0.0895 


-0.0895 



Table 4.5 Static Displacements (in.) 
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From Tables 4.1 - 4.5, it can be seen that the static displacement synthesis 
calculations match those of the Guyan reduction at every step of the process 
on its way to calculate static displacement. 
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V. TIME DOMAIN STRUCTURAL SYNTHESIS 



The following background information on time domain structural 
synthesis is taken from reference [8]. This portion of the thesis will outline 
the methodologies and formulations of the time synthesis technique and the 
classical techniques used to check its accuracy. There will also be a time 
comparison study done in order to quantify the computational efficiency of 
the synthesis technique. 

As mentioned previously, time domain structural synthesis is a 
methodology whereby changes can be made to a given structure, and its new 
transient response determined as the solution of a nonstandard 
nonhomogeneous Volterra integrodifferential equation (VIDE) of the second 
kind. This integral equation is the result of the combination of the impulse 
response functions (IRFs) of the baseline structure, the reaction forces 
generated from the modification to the structure, and the total dynamic 
response of the original system in terms of the convolution integral. This 
greatly differs from the classical method of evaluating a structural change. In 
the classical method, the basic elements which define the structure are 
reconstructed, and then a modal superposition, convolution integral or direct 
integration technique is used on the new structure to solve for the response. 

The motivation for the development of the time domain structural 
synthesis technique stems from the advantages and flexibility that is enjoyed 
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through the use of the frequency domain structural synthesis technique. The 
time domain synthesis can be divided into the same two classes of 
modification and coupling, direct or indirect. Some of the advantages and 
flexibility shared by the two methods are: 1) the formation of an exact 
solution, with the ability to handle damping modifications; 2) arbitrary model 
reduction in the sense that only information about those dofs which are 
needed is retained; 3) the solution phase of the finite element analysis is only 
done once, and the new synthesis model elemental or system matrices are 
never formed. 

Just as with the frequency domain synthesis, this thesis will only focus 
on indirect modifications to a given structure. It will also allow for the 
inclusion of a spring stiffener or a viscous damper between two points, a 
lumped mass on the structure, and implementing a base excitation to the 
structure. A generalized computer program will be presented, which will 
allow for any of the above mentioned changes to be accomplished on some 
baseline structure. The program then returns the transient response for all 
dofs where changes have occurred and any other dofs that may be of interest 
to the user. 

A. GENERALIZED TRANSIENT RESPONSE 

The following equations are for determining the total dynanaic 
response of a system which experiences some sort of excitation. These 
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formulations are taken from references [1 & 2], and are used to verify the 
accuracy of the synthesis after a modification has been made. 



1. Convolution Integral Method 

The convolution integral equation for determining the total dynamic 
response of a system is as follows: 



where: 

{x(t))i, - homogeneous solution which contains the constants of 
integration 

h(t) - matrix of impulse response functions (IRFs) 

F(t) - excitation force 
X - dummy time variable 

The IRF matrix is determined modally using the equation: 



where [O] = the assembled matrix of mass normalized mode shapes. Any 
element of the IRF matrix may be represented as: 



{x(t)} = {x(t)}j, + £[h(t - T)]{F(T)}dr 



( 5 . 1 ) 



[h(t)] = [<D] 



fi)d^t 




( 5 . 2 ) 




G)d t 

r 



( 5 . 3 ) 
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This is a relatively efficient means of calculating a structure's dynamic 
response. However, this method is based on the principal of superposition 
which is valid for linear systems only. Also the modal method of obtaining 
the IRFs does not readily handle non-proportional damping. Therefore, if a 
damping change to the structure is made, or if a non-linear interconnection 
element is used (time domain synthesis with non-linear elements is possible, 
but will not be explored in this study), the response cannot be determined 
using the convolution integral method. 

2. Direct Integration Method 

The direct integration method of calculating a system's dynamic 
response is able to handle the shortcomings experienced by the convolution 
integral method. However, the price is paid in the computational time 
necessary for this solution. For the work illustrated in this thesis, MATLAB's 
ODE45.m fimction is used for this calculation [Ref. 9]. In order to use this 
function, the 2nd order ordinary differential equations (ODEs) of motion: 

[M]{x} + [C]{x} + [KKx} = {F} (3.2) 

must be rewritten as a system of 1st order ODEs. The following is the general 
matrix equation which converts any linear system's equations of motion to a 
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1st order system. 




(5.4) 



where: 



(Zi) = {x} 



{^2} = W 



{zi} = W 

{Z 2 } = {x} 



ODE45.m solves for the vector [{zj { 2 ^]]^ which are the system's response 
displacements and velocities using an adaptive time stepping procedure. 

B. TIME DOMAIN SYNTHESIS FORMULATION 

Consider again the exact same general ndof structure which was 
introduced in Chapter 3 / Section B, and to which structural modifications are 
to be made: 





Figure 5.1 General NDOF Structural System 
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Using the convolution integral, equation (5.1), and the previously 
defined sets, iset, cset, bset, zset, eset for partitioning, the total dynamic 
response of the system can be written as: 



Xi 




Xi 


f* 


1 

1 

-M 


[hicCt-T); 


1 [hj^(t-T)]l' 




>■ = ^ 


Xc 


rl 


0. 

1 


h,^(t-T) 








Xb. 


Jo 

h 


Jhbi(t-'T)] : 


;hb,(t-T); 


1 

1 



Fi(T) 

F,(r) 

F.(t) 



dr 



(5.5) 



Before modification, the force vector refers to externally applied forces 
to the structure. However, following modification, there exists coupling 
forces between the cset and bset dofs. By definition, the iset dofs do not 
experience these coupling forces. Therefore, the partitioned forces may be 
rewritten as: 



(5.6c) 

F, = Ff‘ + F; = Ff ‘ - ( AM,x^ + AC,x; + AK,x; ) (5.6b) 

Fk = Fr + f; = Ff - [4Ct(x; - y)+ AK,(x; - y)] (5.6c) 



where: 

r - the coupling forces due to structure modification. 

[AM], [AC], [AK] - the modification matrices 

X*, x*,x*, - generalized (synthesized) responses after modification. 
Substituting equations (5.6) into equation (5.5), and recognizing that the 
presynthesized responses [Xj x^ x^]^ are those due to the external forces only. 
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the following result is obtained: 



Xil’ 
Xc ^ 



"ju ■- I i^r 
I L I 

J OLlhhlt - l" [hbi,rt - T)J 



[AMJ o' 
0 0 






[ACJ 

0 





(5.7) 



Equation (5.7) is the nonstandard nonhomogeneous VIDE of the second kind. 
In order to obtain the solution for this equation, the integral is discretized 
using a trapezoidal quadrature rule. Also by assuming that there are no 

external excitations to the structure, |F®*‘} = 0, equation (5.7) becomes: 



Xi 

' X, » 

Xb. 



K1.LKI 



J-^bc] ! [-^bb]. 



1 0 




or 



The elemental form of the quadrature matrix is 



o 

1 


0 


0 


1 

o 

o 


f(^). 


f(M 


0 


0 0 


fK). 




f(N)o 


0 0 


• 


* 


• 


0 








- f(N)o. 



(5.8a) 



(5.8b) 



(5.9) 



where the subscript n represents the number of time steps, and the coupling 
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forces equal: 



r[AMj oi|x:i 



(5.10a) 



0 0 0 0 




(5.10b) 




(5.10c) 



Equation (5.8b) can be rearranged and placed in the general form of 
[A]{x*} = {b} where {x*} can be solved for directly using various methods. 
However, due to the possible large size of [A], it is advantageous to leave 
equation (5.8b) as is, and solve for {x*} by iteration. The iteration procedure is 
as follows: 

1) assume the force vector {F^} 

2) calculate the response vector {x^}* 

3) use {xj* to calculate {Xg}*, {Xg}*,and a new {F^} 

4) use the new {FJ and calculate a new {Xg}* 

5) repeat steps 3 & 4 until {Xg}*j,g„ - {Xg}*„jj < convergence tolerance 

Although the time domain structural synthesis is exact in its 
formulation, approximations have now been introduced through the 
solution techniques of the trapezoidal quadrature, and the iteration method. 
Not only is the convergence important, but the stability of this solution is also 
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of concern. The stability and convergence is dependent on the norm of 
and the magnitude of the modifications AM, AC, AK. The norm is defined as 
the largest singular value of the matrix and in this case is driven by the value 
of At. The smaller the value of At, the better the stability and faster 
convergence is reached. The larger the value of the modifications, the greater 
the possibility for instability in the solution. Therefore, due to the solution 
technique, there are now limitations on the modification values. Inherent to 
and dependent on the computer system used, there will be limitations on the 
amoimt of time which can be covered by this analysis due to memory 
capabilities. 

C TIME DOMAIN SYNTHESIS COMPUTER CODE 

The goal in the formulation of the time domain synthesis computer 
programs is to perform the synthesis for an indirect modification operation to 
a given structure, allowing the same type of modifications presented 
previously in the frequency domain synthesis chapter of the thesis. The 
program then returns the transient response for all dofs where changes have 
occurred and any other dofs that may be of interest to the user. Appendix C 
shows the computer codes used to perform the time domain synthesis, and 
perform a comparative analysis of the synthesis versus classical methods. 

The programs presented in Appendix C can be segregated into two 
categories. The first is the calculation of the synthesized transient responses 
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due to a base excitation on the structure, and their comparison to the 
classically calculated responses. One of the programs in this category uses the 
convolution integral method in order to check the synthesis, while the other 
uses direct integration. The second category is the calculation of the 
synthesized transient responses where the different aspects of the synthesis 
technique have been modularized into MATLAB functions. The following is 
a diagram of the flowpath of the modularized time domain synthesis 
programs: 




Figure 5.2 Flowpath of the Modularized Time Synthesis Programs 



A more detailed description of each function is contained within Appendix C. 
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D. ILLUSTRATIVE EXAMPLES 



The purpose of this section is to illustrate the use of the time domain 
synthesis technique and validate its accuracy. A computational time 
comparison will be done between the synthesis and classical methods in order 
to assess the computational efficiency of this technique. The three structural 
systems previously described, mass-spring-damper, beam and mass-plate, are 
used to accomplish this goal. 

1. Mass-Spring-Damper System 

Consider the following mass-spring-damper system: 



Figure 5.3 Mass-Spring-Damper System Experiencing Base Excitation 

The solid elements represent the original structure and the dashed elements 
represent modifications. Although not shown in the figure, the original 
structure is proportionally damped using [C] = a[K]. The modifications shown 
are arbitrarily selected. 




k = 1000 Ib/in 



m = 0.259 Ib-sVin 



Ak 24 = 250 IbAn 
Akb = 2501b/in 




m 
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a. Spring Modification Only 

The following figure is the resultant synthesized and classically 
determined transient response XjCt) which represents the response at dof 2 
due to the base excitation. For this example neither the mass nor the damper 
modification was included. 



Transient Time Response at dof 2 




Figure 5.4 Transient Response Mass-Spring-Damper System 

From the graph produced, it can been seen that the responses are identical, 
proving that the synthesis technique is exact, and that the computer coding 
for the arbitrary changes is correct. These results were produced using the 
tsynconv.m program. 
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bt Spring & Mass Modifications 

The previous example is now repeated but with the mass change 
at dof 2, Amj = .02 Ib-s^/ in included. The following figure shows the results of 
this additional modification. 

Transient Time Response at dof 2 




Figure 5.5 Transient Response Mass-Spring-Damper System w/Spring-Mass 

Modification 

From the graph produced, it can been seen that the mass modification had no 
effect on the accvuacy of the solution. However, the addition of the mass 
modification did have a large impact on the computational time required for 
the s)mthesis method. This will be discussed in more detail later. 
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c Spring & Damper Modifications 

The spring modification only example is now repeated but with 
a damper change at the base Ac^ = 2 Ib-s/in included. The following figure 
shows the results of this additional modification. 



Transient Time Response at dof 2 




Figxire 5.6 Transient Response Mass-Spring-Damper System w/Spring- 

Damper Modification 

From the graph produced, it can been seen that the responses are identical 
with respect to the resolution of the plot graphics. Also, to further validate 
the correctness of the computer algorithms, the effect of the added damper 
can be seen when comparing Fig. 5.4 & Fig. 5.6. Since the non-proportional 
damper change was made to the structure, the direct integration method was 
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used to obtain the classical solution. Although there was no real effect in the 
accuracy of the solution, the effect in computational time is profound for the 
classical method. These results were produced using the ts)mode45.m 
program. 



d. Spring , Mass, & Damper Modifications 
The example is now repeated with all modifications included. 
The following figure shows the results of these modification. 



Transient Time Response at dot 2 




Figure 5.7 Transient Response Mass-Spring-Damper System w/ Spring-Mass- 

Damper Modification 
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Just as in all the other cases, it can been seen that the responses are extremely 
close to identical. These results were also produced using the tsynode45.m 
program. 

2. Beam System 

Now consider the following beam system: 



Akb = 1000 Ib/in 
Acijo = 2 Ib-s/in 
y(t) = Yo(e“^ - 
Yo=l 



Figure 5.8 Beam System Experiencing Base Excitation 

The beam parameters and finite element modeling of the above structure is 
identical to that used in the frequency domain synthesis chapter of the thesis. 
The modifications shown are arbitrarily selected. 

a. Base Spring Modification Only 

The following figure is the resultant synthesized and classically 
determined transient response x„(t) which represents the response at dof 11 
due to the base excitations. As can be seen from Fig. 5.7, this is not a dof where 
change has occurred or where the base excitation is attached, but just a dof 




✓ X ^ 
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where there is interest in the response. For this example the damper 
modification Ac^jo was not included. 



Transient Time Response at dot 11 




Figure 5.9 Transient Response Beam System 



From the graph produced, it can been seen that the responses are not 
identical, but follow each other with relatively little error. The difference 
between the two methods can be attributed to the numerical solution of the 
synthesis method. As the time step for analysis is decreased, the numerical 
solution approaches the exact solution. However, the use of smaller time 
steps is restricted by computer memory limitations. 
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b. Damper Modification 

The previous example is now repeated but with the damper 
change at dof 13, = 2 Ib-s/in included. The following figure shows the 

results of this additional modification. 

Transient Time Response at dof 1 1 




Figure 5.10 Transient Response Beam System w/Damper Modification 

From the graph produced, it can been seen that the responses are not 
identical, but follow each other with relatively little error. The same line of 
reasoning as in the undamped case is applicable here also. A smaller time step 
yields more accurate data but, memory limitations inhibits its usefullness. 
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3. 



Mass-Plate System 



Now consider the following mass-plate system: 




Cb y(t) ,Xs 



Figure 5.11 Plate-Mass System Experiencing Base Excitation 

The plate parameters and finite element modeling of the above structure is 
identical to that used in the frequency domain synthesis chapter of the thesis. 
The modifications shown are arbitrarily selected, but are in accordance with 
the changes which will be made during an optimization problem. 

a. Spring Modifications Only 

The following figure is the resultant synthesized and classically 
determined transient response Xj^^Ct) which represents the response at dof 109, 
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the computer, due to the base excitations at the plate's comers. For this 
example none of the damper modifications are included. 



Transient Time Response at dof 109 




Figure 5.12 Transient Response Mass-Plate System 



From the graph produced, it can been seen that the responses are not identical 
but follow each other with relatively little error. The error difference between 
the two methods can be attributed to the numerical solution technique, and 
the step size of the time vector. As the step size decreases, the synthesis 
solution approaches the exact solution. However, memory limitations on the 
computer system used for this research, prohibit meaningful use of a smaller 
time step. 
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bi Spring & Damper Modifications 

The spring modifications only example is now repeated but with 
the prescribed damper changes at the base and between the computer and the 
plate. The following figure shows the results of these additional 
modifications. 



Synthesized Transient Time Response at dof 1 09 




Figiure 5.13 Transient Response Mass-Plate System w/ Spring-Damper 

Modifications 

Due to the size of this model and the damper modifications being made, it 
was not possible to compare the synthesis solution with a classical solution. 
The computational time required for a classical solution would be great. 
However, from comparing Figs. 5.11 & 5.12, the slight effect of the additional 
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damping can be seen in the peak amplitude of the response, giving some 
credence to the solution. Also, the other previously run damper modification 
examples show that this technique does in fact work. 



4. Computational Time Comparisons 

The following table is a compilation of all the computation times that 
it took for the previously presented examples to complete its process. 







Synthesis 


Classical Methods 








(secs) 


Convolution 

(secs) 


Integration 

(secs) 


At 

(secs) 


Mass- 

Spring 


No Damper 
No Mass 


15.62 


1.28 




.001 




Damper 
No Mass 


23.12 




379.32 


.001 




Mass 

No Damper 


194.05 


1.27 




.001 




Mass 

Damper 


197.89 




339.55 


.001 


Beam 


No Damper 


107.38 


8.56 




.0003 




Damper 


183.75 




5687.88 


.0003 


Mass-Plate 


No Damper 


27.39 


15.77 




.0008 




Damper 


99.99 




t — > oo 


.0008 



Table 5.1 Synthesis vs. Classical Computational Times 
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From the table, it can be observed that for smaller structures, the classical 
method seems more efficient. This is due to the fact that the advantage of the 
iteration method is lost when dealing with small structures. The 
computational time for the iteration solution, or in other words the synthesis 
technique, is driven by the magnitude of the time step and the number of 
time steps. On the other hand, the classical methods' solution times are 
driven by the model size. Therefore, depending on the time step size, a very 
small and very simple model may take a very long time to solve using time 
domain synthesis. 

Another interesting aspect of the comparison is how sensitive the 
synthesis technique is to mass changes. Since the mass change does not add 
additional dofs to the model, the classical method's computational time is not 
increased by this modification. On the other hand, a mass change for the 
synthesis method constitutes the use of the calculated second derivative of 
the synthesized response. This derivative is accomplished using a finite 
differencing scheme and inherently adds the possibility of additional 
instability to the solution. 

As previously alluded to, the finite differencing could be the reason 
why the synthesis method is greatly affected by mass changes. However, 
damper changes, which also prompt the use of finite differencing, did not 
substantially increase computational time for the synthesis method. This is in 
contradiction to the thought that the increased instability is caused by the 
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finite differencing scheme. Another factor which plays into the convergence 
and stability of the solution is the magnitude of the modification. Relative to 
each other and this problem, the mass change may be of greater magnitude 
than the damper change and therefore costs more in computational time. 

The largest and most noticeable disparity is in the time required to 
classically calculate the dynamic responses following a damper modification. 
The direct integration method is very inefficient and becomes basically 
unusable for very large complex models, or in the use of an optimization 
routine. 
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VI. OPTIMIZATION 



The following background information and the motivation for this 
chapter can be attributed in part to reference [7]. Inherent to the operation of 
military vehicles, whether they operate on land, sea or in the air, is the 
presence of a shock and vibration environment. Computer and electronic 
equipment located on these vehicles can malfunction or be damaged as a 
result of the severe effects of this environment. Therefore, the need arises for 
a properly designed isolation system, which will absorb the shock and 
vibrations and therefore protect the equipment. This is an especially 
important concept in today's military with the relatively recent practice of 
using Commercial-Off-The-Shelf (COTS) equipment vice the traditionally 
expensive shock and vibration hardened military equipment. 

In general, it is relatively straightforward to design an isolation system 
to protect against shock only or vibration only. However, these two designs 
are in direct conflict with each other [Ref. 1]. Therefore, in the presence of 
both types of excitations and various other constraints, the optimal design of 
an isolation system becomes more complicated. The need then arises for 
some sort of optimization scheme which will minimize the response of the 
component to be protected, while simultaneously satisfying the constraints 
which have been placed on the system. 
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The complexity of this problem does not necessarily end here. If the 
structure to be isolated and the isolation system to be designed is large and 
complex, the analysis process to calculate the desired responses can be very 
computationally demanding. Furthermore, if there are multiple design 
variables which may be altered in order to achieve the desired result, the 
search patterns for the optimal solution are more complex, which leads to 
more iterations of the optimization process. The combination of these two 
factors can lead to the conclusion that the optimal design of a large complex 
isolation system is impractical. 

With the previously presented arguments in mind, the use of 
computationally efficient synthesis techniques for the calculation of the 
system's response is the obvious answer. Since the structure is large and 
complex, analysis time could be relatively long. The larger the system, the 
greater the number of possible design variables. This means that the search 
for the optimal solution is even more complex, and will require more 
iterations. This portion of the thesis will demonstrate the use of the 
previously presented structural synthesis techniques in the optimal design of 
a shock and vibration isolation system. A general computer algorithm will be 
formulated and used in order to illustrate the use of structural synthesis in 
optimization design. 



82 



A. 



GENERAL OPTIMIZATION PROBLEM FORMULATION 



In an optimization problem, the quantity to be maximized or 
minimized is termed the objective function. The quantities which may be 
altered in order to find this optimal value are termed the design variables. 

The constraint fimctions are also a function of the design variables, and three 
different types of constraints may exist in an optimization problem. The first 
type is an inequality constraint, where the value calculated from the 
constraint function is less than or equal to some prescribed limit. The second 
type is an equality constraint. Here the value calculated from the constraint 
function must be equal to some prescribed value. The third constraint type is 
side constraints. This is when upper and lower limits are placed on the design 
variables. In a constrained optimization problem, the objective fimction is 
maxiiiuzed or minimized, while ensuring that the user defined values of the 
constraint functions are not violated. 

For the design of an isolation system, the minimization of the 
equipment's response over some frequency range, due to vibration, could be 
the main concern. At the same time, the entire structure's response due to a 
shock input, and the static displacement or 'sag' of the system could be 
constrained. The physical parameters which can be modified to achieve this 
goal could be the stiffness and damping values of the isolators. The design of 
this isolation system can be converted to the following optimization problem: 
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Minimize (Objective Function): 

Computer response due to random vibration 
Subject to (Constraints): 

Maximum static displacements < limit 
Maximum dynamic isolator stretch due to shock < limit 
Maximum computer acceleration due to shock < limit 
Min and max changes in the isolator stiffness values < limit 
Min and max changes in the isolator damping values < limit 
Design Variables: 

Changes in the isolator stiffness values 
Changes in the isolator damping values 

1. Calculation of Objective Function 

Since the objective fimction is based on the response due to random 
vibrations, the input base excitation is in the form of a power spectral density 
(PSD) fimction. The response is therefore calculated as the mean square value 
in terms of the system response function (FRF), and the PSD of the input. The 
equation for the mean square value of the response is given as [Ref. 10]: 



From equation (6.1), it can be seen that the most accurate means of 
minimizing this equation would be if the input PSD, Sif), is known. Since the 




( 6 . 1 ) 
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goal is to formulate a general optimization program, it is not inconceivable 
and is appropriate to assume S(/) to be constant over some frequency, /, < / < 
/j. Therefore, minimizing the mean square response over the appropriate 
frequency range is equivalent to minimizing: 



It can be seen that for every iteration of the optimization process, the 
calculations involved in equation (6.2a) must be performed. If the FRFs of the 
baseline structure is provided or calculated prior to the iterations begin, the 
frequency domain structural synthesis technique may be used to calculate all 
subsequent FRFs and the objective function is properly written as: 



2. Calculation of Constraints 

The first constraints mentioned are the maximum static displacements 
of the structure. Since these values must also be recalculated every iteration, 
an efficient means of doing so needed. Various reasons make this problem a 
prime candidate for the use of static displacement synthesis: 

• the static displacement at select locations vice the entire structure is 
desired, therefore calling for some type of reduction or 
condensation technique. 

• the elemental and system matrices of the structure may not be 
available. 




(6.2a) 




(6.2b) 
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• the modal solution of the eigenvalue problem has already been 
accomplished. 

• the isolator stiffness, which are design variables, act in vertical 
translation direction only (isolator damping has no effect on static 
problem). 

Another constraint is the maximum dynamic isolator stretch due to 
the shock input. In order to calculate this constraint, the transient response of 
the structure at the dofs where the isolators are connected must be obtained. 
The displacement of the base input as a function of time must also be known. 
Since the change in the isolator damper represents nonproportional 
damping, the direct integration method would be used to solve this problem. 
This calculation would be terrible, in reference to computing time for a single 
iteration of the optimization process. Despite other reasons, this alone is 
enough motivation for the implementation of the time domain synthesis 
technique for the calculation of the dynamic responses. The last constraint, 
the computer acceleration, is obtained by simply applying a finite differencing 
scheme to the computer's displacement time history. All the above 
mentioned constraints are known as inequality constraints, and will satisfy 
the general normalized form gj(X) < 0 [Ref. 11]. 
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3. Design Variables 

There are limits on the amoimt in which the isolator stiffness and 
damping can be changed. For example, the lower bound for the change in an 
isolator stiffness should not be equal to the negative of the original value of 
the isolator. If it is, the optimization program might use this value and this 
constitutes the total removal of the isolator and changing the basic structure 
of the system. The upper bound for the change is limited by what is physically 
realistic. These type of constraints on the design variables are known as side 
constraints [Ref. 11]. 

The relative value of design variables to one another is very important 
when considering the efficiency and reliability of the numerical optimization 
process [Ref. 11]. The function space of the problem, which is defined by the 
design variables, should be as symmetric as possible. When looking at the 
design variables, isolator stiffness and isolator damping, it can be seen that 
effective changes to isolator stiffiiess are on the magnitude of thousands, 
while the isolator damper changes are on the magnitude of ones. This large 
disparity creates a very nonsymmetric function space. In order to alleviate 
these problems, the design variables are scaled so that their magnitudes are 
comparable. A symmetric function domain allows for the faster 
determination of the optimization by decreasing the number of iterations 
necessary to find the optimal solution[Ref. 11]. 
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B. OPTIMIZATION COMPUTER CODE 



The optimization tool used to solve this problem is the MATLAB 
optimization toolbox l.Od. More specifically, the function constr.m, which is 
designed for optimization of a constrained, nonlinear, multivariable 
problem, is used. This function uses the Sequential Quadratic Programming 
(SQP) method in order to solve this problem [Ref. 12]. The SQP method is 
used in order to determine search directions for the design variables at every 
iteration [Refs. 11 & 12]. One of the limitations of this, and other traditional 
nonlinear optimization codes, is that they may only give local solutions. 
Also, when the problem is infeasible, constr.m attempts to minimize the 
maximum constraint value. What this can lead to is, the problem iterating to 
the user defined iteration limit, the search pattern extending outside the 
bounds of the design variables, and no optimum solution ever being found. 
Therefore, the user must have good understanding of the problem and the 
constraints imposed upon it. 

The goal in the formulation of the optimization computer programs is 
to provide a means of determining the optimal changes in isolator values, in 
order to minimize some dynamic response, while simultaneously satisfying 
all prescribed constraints. These programs first allow for the loading of a 
general finite element model of the structure. This baseline structure is given 
in terms of its FRFs and IRFs, or in terms of its system matrices from which 
the baseline FRFs and Ifs must be calculated. Modifications to the baseline 



88 



struchire which will act as design variables are then designated. User defined 
initial values of design variables and constraint values are designated. During 
the optimization process, all structure responses will be calculated using 
structural synthesis techniques. With minor coding adjustments, the ability 
to interchange objective function and constraints exists. For example the 
optimization program could be the minimization of the computer's 
acceleration due to shock, while the response to a random vibration input is a 
constraint. Appendix D shows the computer codes used to perform this 
process. The following figure is a diagram of the flowpath of the optimization 
program. 
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Figure 6.1 Flowpath of the Optimization Program 
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C ILLUSTRATIVE EXAMPLE 



The purpose of this section is to illustrate the use of the structural 
synthesis techniques in the optimization of a shock and vibration isolation 
system. The figure below shows the mass-plate structural system which was 
used in the previous examples and is also used here. 




1. Problem Formulation 

The following is the formal problem statement of the optimization 
problem: 

Minimize: F(Aki,,Ak^,ACb,Ac^) = £^|H*(/)| d/ 
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Zp_stat 



- 1<0 



Subject to: 



gi = 



g2 = 



g3 = 
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Ss “ 



- 1<0 



max Zp_stat 

z^_stat 
max z^_stat 



kp_stretch 
max kp_stretch 

k^,_stretch 
max kj_stretch 



Z£_accel 
max z^_accel 



- 1<0 



— 1<0 



- 1<0 



^000 <Akb <5000 
^000 <Ak,< 5000 



0 < ACb < 5 
0 < Ac, < 5 



where: 

H*(/) = frequency response function of the computer 
Zp_stat = plate static displacement [in] 

z,_stat = relative static displacement between computer and plate [in] 
maxZp_stat = maximum allowable plate static displacement [in] 
maxz,_stat = maximum allowable static displacement between 
computer and plate [in] 

kp_stretch = isolator stretch between base and plate due to shock [in] 
k^stretch = isolator stretch between plate and computer due to shock 
[in] 

max kp_stretch = maximum allowable isolator stretch between base and 
plate [in] 

max k,_stretch = max allowable isolator stretch between plate and 
computer [in] 

z^accel = absolute acceleration of computer due to shock [in/s^] 
max z,_accel = max allowable absolute acceleration of computer [in/s^] 
Akb, A k„ ACb, Ac, = the design variables. 
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The following is information that must be provided to the 
optimization program for its execution, and for the imderstanding of the 
user. The following table lists the initial, lower boimd, and upper bound 
values for the design variables. 





Akb 


Ak, 


ACb 


Ac, 


initial value 


1000 


1000 


1 


1 


lower bound 


-4000 


-4000 


0 


0 


upper bound 


5000 


5000 


5 


5 



Table 6.1 Optimization Inputs 



From the table it can be noticed that an initial value for the design variables is 
used. Since the optimization routine locates local minimum, the start 
position of the search is very important. The baseline values for the isolator 
elements which are to be modified are: 



initial = 10000 Ib/in 

initial = 5000 Ib/in 

initial c^ =0 Ib-s/in 
initial c^ =0 Ib-s/in 



The constraint values for the inequality constraint are: 

maxz^stat = -0.08 in 

p- 

maxz^stat = -0.03 in 
max kp_stretch = 1.0 in 
max k^_stretch = 1.0 in 
max z^accel = 30g's 
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2. Optimization Results 

The previously described problem is executed using the isomain.m and 
isosub.m programs. During the execution of the optimization program, 
output information is provided through the MATLAB Command Window 
(Table 6.2). 



f-COUNT 


FUNCTION 


MAXfg] 


STEP 


Procedures 


5 


36012.1 


1.40603 


1 




10 


15494 


0.0364448 


1 




15 


6969.56 


-0.0262839 


1 




22 


6691.27 


-0.0370827 


0.25 




29 


6372.67 


-0.0682749 


0.25 




34 


6372.59 


-0.0693421 


1 




39 


6370.5 


-0.0631356 


1 


mod Hess(2) 


44 


6370.49 


-0.0687934 


1 


mod Hess 


49 


6370.49 


-0.0689084 


1 


mod Hess(2) 


54 


6370.48 


-0.067107 


1 




59 


6370.47 


-0.0527099 


1 


mod Hess 


64 


6370.45 


-0.0667607 


1 




69 


6370.44 


-0.0687255 


1 


mod Hess(2) 


74 


6368.61 


-0.0685718 


1 


mod Hess(2) 


79 


6368.61 


-0.0688113 


1 


mod Hess 


80 


6368.61 


-0.0667838 


1 


mod Hess 



Table 6.2 Optimization l.Od Generated Output 



Once the optimization program terminates successfully, post 
processing occurs and the following figures and summary information is 
provided. The first figure is a graph of the value of the design variables at 
each iteration. 



94 



Change in Variables vs. iterations 




Figvire 6.3 Values of the Design Variables 



The next figtire is a graph of the absolute value of the isolators at each 
iteration of the optimization process. 
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Change in Isolation Constants vs. Iterations 
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Figure 6.4 Absolute Values of the Isolator Elements 



The next figure is a graph of the FRF of the computer before the optimization 
begin, and at its termination. 
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Figiire 6.5 Computer FRF Before and After Optimization 



The effect of the changes to the structure, particularly the damping additions, 
can definitely be seen in this figure. The objective function is defined as the 
area under this curve after it is squared. Therefore, by lowering the peak 
values of the response, the objective function is ultimately lowered. The next 
figure is a plot of how the value for the objective fimction changed as a 
function of the optimization iterations. 
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Figure 6.6 Values of the Objective Function 

From Fig. 6.6, the dramatic effects of the optimization can really be seen. The 
initial value of the objective function was 36012.1 and after the isolator 
modifications the objective function equals 6368.31. The following 
information is a summary of the optimization process and is displayed at the 
MATLAB Command Window: 

Optimization Terminated Successfully 

Active Constraints: 
ans = 3 

The recommended change in base isolator stiffness (lb /in) is — > 
final_del_kb = 2.3271e+03 Ib/in 
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The recommended change in computer isolator stiffness (lb/ in) is — > 
final_del_kc = 4.0849e+03 Ib/in 

The recommended change in base isolator damping (Ib/in/s) is — > 
final_del_cb = 2.7678 Ib-s/in 

The recommended change in computer isolator damping (Ib/in/s) is — > 
final_del_cc = 4.6266 Ib-s/in 

Constraints: 

gi = -0.1814 g2 = -0.6334 ga = -0.5237 g4 = -0.6981 gg = -0.0668 
elapsed_time = 9.5550e+03 secs = 2hrs 39min 

It can be seen from the presented information that it took the 
optimization program 80 iterations to recommend the presented changes to 
the structure in order to optimize its response against random vibrations, and 
meet the established constraint criteria. The most notable portion of this 
example is that it accomplished this task on a 109 dof structure, which 
underwent non-proportional damping changes at every iteration, in only 2 
hours and 39 minutes. In other words, this program performed 80 transient 
response analyses, calculated the static displacement 80 times, and calculated 
the frequency response fimction and the integral of that curve 80 times, for a 
109 dof structiire. 
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VII. CONCLUSIONS / RECOMMENDATIONS 



The main objective of this thesis was to investigate the use of 
structural synthesis techniques as rapid reanalysis procedures in the optimal 
design of large shock and vibration isolation systems. In order to use these 
synthesis techniques liberally, an investigation into the accuracy of the 
methods, and the computer algorithms formulated, had to be accomplished. 

It is important to mention at this point that the time and frequency structural 
s)mthesis methods used in this study, primarily covers base excitations. This 
is a new extension of the more familiar synthesis formulations. The static 
displacement synthesis, is however a completely new development. These 
formulations had to be general enough to handle different structures and to 
handle various types of modifications to the structure. The synthesis 
techniques were then coupled with an optimization routine, and together 
they were used in finding the optimal changes to a given isolation system. 

A. CONCLUSIONS 

From the various analyses conducted, the following conclusions can be 
drawn: 

1. The frequency domain structural synthesis technique is a very 
efficient and accurate method of calculating the new FRFs of a 
structure following indirect modifications. It is nonrestrictive and 
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arbitrary as to where and what type of modifications can be made to 
the structure. 

2. The static displacement structural synthesis is also a very efficient 
and accurate method of calculating the new static displacement of a 
structure following modifications. The accuracy of this method does 
depend on the use of a finite differencing scheme which in Uurn is 
dependent on the variable step size. It is also found that this method 
is restrictive in its use. Certain criteria, which are identified in 
Chapter IV of this work, must be met in order that this synthesis 
technique can be applied properly. For the isolation problems 
presented in this work, the required criteria were met. 

3. The time domain structural synthesis technique is an efficient and 
accurate method of calculating the new transient response of a 
structure following indirect modifications. It is nonrestrictive and 
arbitrary as to where and what type of modifications can be made to 
the structure. However, since numerical techniques are used to 
solve for the responses, some error of approximation is entered into 
the solution. As the time step is decreased, the numerical solution 
approaches the exact solution. However, current solution 
algorithms are restricted by memory requirements. 

4. In finding the optimal design of a relatively large shock and 
vibration isolation system, the use of the synthesis techniques 
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proved to be an invaluable asset in the optimization routine. With 
the use of the synthesis techniques, in 2 hours and 39 minutes, it 
was possible to perform 80 FRF, static displacement, and transient 
analyses on a 109 dof system, which was non-proportionally 
damped. To get an idea of how efficient this is, recall fi-om the time 
domain synthesis analysis (Chapter V), that it took approximately 1 
and 1/2 hours to calculate one transient analysis of a 22 dof beam 
structure. 

B. RECOMMENDATIONS 

Throughout this study, some weaknesses of the techniques, and 
opportunities to enhance the techniques were touched upon. 
Recommendations to address these weaknesses, and take advantage of these 
opportunities include: 

• All modal analyses were done to include all modes of the structure. 
A modal truncation criteria should be established, and studies 
should be done to ensure that these criteria hold up for synthesis 
also. This way it would be possible to take advantage of time saving 
aspects which are inherent to modal superposition techniques. 

• Use actual FRF and IRF data from a real structure, and predict 
structural responses due to modifications. 
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• Build CAD models to validate damping modifications to large 
structures. 

• Investigate to obtain a general rule for choosing the optimum AQ to 
be used in finite difference scheme to calculate static response. 
Another option would be to investigate the use of the closed form 

2nd derivative of H* to obtain ^ . 

• Implement relaxation schemes for the iterative solution in time 
domain synthesis. Also a study should be done in order to 
determine what the proper time step size should be to ensure 
convergence and stability. 
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APPENDIX A. FREQUENCY DOMAIN SYNTHESIS COMPUTER 



CODES 



The following is a list and a brief description of the main MATLAB 

computer codes that were written in order to perform the frequency domain 

synthesis calculations. 

• fsynstr.m - performs the frequency domain synthesis on a structure which 
experiences an external force excitation directly on the structure, and 
compares this solution to the FRF calculated using classical methods. 

• fsynbase.m - performs the frequency domain synthesis on a structure 
which experiences a base excitation, and compares this solution to the FRF 
calculated using classical methods. 

• fsynmain.m - the main program which calls the modularized frequency 
synthesis programs and does the post processing of the results. 

• change.m - loads the baseline structure to be modified, allows for 
modifications to the baseline structure, and calls the fimction delta.m 

• delta.m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

• modal.m - solves for the structure's natural frequencies, mode shapes and 
other modal matrices and vectors. 

• frfmodal.m - calculates the structure's FRFs. 
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• frfsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized FRFs. 

• frf sift.m - selects the FRF that the user wants to perform post processing 



on. 

The full codes are contained on subsequent pages. 
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fsynstr.m 



% This program is designed to calculate FRFs by inverting the impedance matrix. 

% Changing the mass, stiffness, and damping matrices, new FRFs will be calculated by 
% reformulating the matrices, and using Ae synthesis method in the frequency domain. 
% 

clear 

% Formation of the elemental matrices 
% 

nel=4; 

nn=nel+l; 

dof=nel+l; 

ndk=l; 
knel(l)=0; 
for spring=l:ndk; 
kcon(spring)=100; 
knel(spring+l )=4; 

for ele=knel(spring)+l :knel(spring+l) 
kcon_ele(ele)=kcon(spring); 
end 
end 

ndm=l; 

mnel(l)=0; 

formass=l;ndm; 

m(mass)=l 00/386.4; 
mnel(mass+l)=5; 

for ele=mnel(mass)+l :mnel(mass+l) 
mele(ele)=m(mass); 
end 

end 

% Formulation of global stiffness, damping and mass matrices from elemental matrices 
% 

Kglo=zeros(dof,dof); 

Cglo=zeros(dof,dof); 

Mglo=zeros(dof,dof) ; 
for n=l:nel 

kele=kcon_ele(n)*[ 1-1 
-1 1 ]; 

rcl=n; 

rc2=n+l; 

Kglo(rcl;rc2u-cl;rc2)=Kglo(rcl:rc2,rcl:rc2) + kele; 
end 

Mglo=diag(mele); 
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alpha=0; beta=.0005; 
Cglo=alpha*Mglo + beta*Kglo; 



bc=2; 

ifbc=l 

K=Kglo(2:dof,2:dof); 

K((dof-l),:)=[]; 

K(:,(dof-l))=[]; 

M=Mglo(2:dof,2:dof); 

M((dof-l), :)=[]; 

M(:,(dof-l))=[]; 

C=Cglo(2:dof,2:dof); 

C((dof-l),:)=[]; 

C(:,(dof-l))=[]; 
ndof=dof-2; 
elseif bc=2 
K=Kglo(2:dof,2:dof); 

M=Mglo(2:dof,2:dof); 

C=Cglo(2:dof,2:dof); 
ndof^of-1; 
elseif bc=3 
K=Kglo(2:dof,2:dof); 

K((dof-l), :)=[]; 

K(:,(dof-l))=[]; 

M=Mglo(2:dof,2:dof); 

M((dof-l),:)=[]; 

M(:,(dof-l))=[]; 

C=Cglo(2:dof,2:dof); 

C((dof-l),:)=[]; 

C(:,(dof-l))=[]; 

ndof=dof-2; 

else 

K=Kglo; M=Mglo; C=Cglo; 

ndof:^iof; 
end 



% This portion of the program is designed to accept changes to the original 
% structure and recalculate the FRF using classic^ metiw^ of reforming [M], 
% [K], and [C] matrices using the inverse of the impedence matrix. Also 
% recalculates the FRF using structural synthesis in the frequency 
% domain. 

% 

change_model=l; 
while change_model=l; 

K_c=K; M_c=M; C_c=C; % initializes change matrices to 

% original matrices 

changek=2; changec=2; changem=2; % initializes logic variables 

% to 'no' status 
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chk=0; chc=0; chm=0; % initializes counter for # of changes 

% to matrices 

kdofl=[l; kdof2=Q; % initalizes the matrices which keep 

cdofl=[]; cdof2=0; % track of the dofs where changes occur 

mdof=D; 

cset=0; % initializes cset matrix to zero to avoid null index error in 

% handling the spring to ground if no changes are made to the 
% structure. 

changek=inputCWould you like to change your stiffness matrix? l=yes 2=no ')> 
while changek—1 

chk=chk+l; % counter for determining # changes to stiffness matrix. 
lk(chk)=input('input value of added spring (Ibs/in) ’); 

kdofl(chk)=input(’input constrained dof where 1st end of added spring is applied '); 
kdof2(chk)=input('input dof where 2nd end of added spring is applied or 0 for ground 

’); 

% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if kdof 1 (chk)~=cset 
cset=[cset kdofl(chk)]; 
end 

if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 

% Classical method of reformulatting the [K] matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk)) -h Ik(chk); 
if kdof2(chk)~^ 

K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)dcdof2(chk)) - lk(chk); 
K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)dcdofl(chk)) - Ik(chk); 
K_c(kdof2(chk)dcdof2(chk))=K_c(kdof2(chk)dcdof2(chk)) + Ik(chk); 
end 

changek=input('Are there any other changes to the stiffness matrix? l=yes 2=no '); 
end 

changec=inputCWould you like to change your damping matrix? l=yes 2=no '); 
while changec==l 

chc=chc+l ; % counter for determining # changes to damping matrix. 

flag_c=l; % mapping matrix flag 
lc(chc)=input(’input value of added damper (Ibs-sec/in) ’); 

cdofl(chc)=input('input constrained dof where 1st end of added damper is applied '); 
cdof2(chc)=input(’input dof where 2nd end of added damper is applied or 0 for 
ground ’); 
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% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if cdofl(chc)~=cset 
cset=[cset cdofl(chc)]; 
end 

if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 



% Classical method of reformulatting the [C] matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc)) + lc(chc); 
if cdof2(chc)~^ 

C_c(cdofl(chc),cdof2(chc))=C_c(cdofl(chc),cdof2(chc)) - lc(chc); 
C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc)) - lc(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 

changec=input(’Are there any other changes to the damping matrix? l=yes 2=no '); 
end 

changem=inputCWould you like to change your mass matrix? l=yes 2=no ’)» 
while changem=l 

chm=chm+l ; % counter for determining # changes to mass matrix. 

lm(chm)=input(’input value of added mass (lbs-sec'^2/in) '); 
mdof(chm)=input(’input constrained dof where mass is added ’); 

% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 

% Classical method of reformulatting the [M] matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + hn(chm); 

changem=inputCAre there any other changes to the mass matrix? l=yes 2=no *); 
end 

% Handling of spring added to ground 
% 

ground = find(cset=0); 
cset(ground) = D; 
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% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 

% 

iset=inputCWhat dofs other than those where change occured, are you interested in? '); 

% Formation of eset matrix which is i U c. 

% 

eset=[iset cset]; 

% Choosing the FRF interested in 
% 

resp_dof=inputCWhat response dof are you interested in? '); 
inp_dof=input('What input dof are you interested in? '); 
frf_index=find(eset=resp_dof I eset==dnp_dof); 
if resp_dof==inp_dof 
fif_index=[frf_index fif_index]; 
end 

% Structural Synthesis method of calculating FRF 



% Calculation of unchanged FRF 
% 

freq=.01:.01:10; 
omega=2*pi*freq; 
j=length(omega); 
for w=l:j 

Z=K+i*omega(w)*C-omega(w)'^2*M; 

H=inv(Z); 



% Rearrangement of original [H] to form [Hee], [Hec], [Hcc], [Hce] 
% Also provides for if user makes no changes to the model 
% 

Hee=H(eset,eset); 

if cset=n 

Hec=0; Hcc=0; Hce=0; 

else 

Hec=H(eset,cset); 

Hcc=H(cset,cset) ; 

Hce=H(cset,eset); 

end 

% Formation of [del_K], [del_C], [del_M], & [del_Z] 

% 

del_K=zeros(length(cset)); % intializes the change matrices 
del_C=zeros(length(cset)); % to zero and sets their sizes 
del_M=zerosOength(cset)); % as [cset x cset] matrix 

if chk~=0 
for a=l:chk 

indexkl=find(kdofl(a)==cset); % finds where in cset matrix 
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indexk2=find(kdoG(a)=cset); % the change in k occurs 

del_K(indexkl,indexkl)=del_K(indexkl4ndexkl) + lk(a); 
del_K(indexkl,indexk2)=del_K(indexkl,indexk2) - lk(a); 
del_K(indexk2,indexkl)=del_K(indexk2,indexkl) - lk(a); 
del_K(indexk2,indexk2)=del_K(indexk2,indexk2) + lk(a); 
end 
end 

if chc~=0 
fora=l:chc 

indexcl=find(cdofl(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); % the change in k occurs 

del_C(indexcl,indexcl)=del_C(indexcl4ndexcl) + lc(a); 
del_C(indexcl,indexc2)=del_C(indexcl,indexc2) - lc(a); 
del_C(indexc2,indexcl)=del_C(indexc24ndexcl) - lc(a); 
del_C(indexc2,indexc2)=del_C(indexc2,indexc2) + lc(a); 
end 
end 

if chm~=0 
for a=l:chm 

indexm=find(mdof(a)==cset); % finds where in cset matrix 

% the change in k occurs 

del_M(indexm,indexm)=del_M(indexm4ndexm) + lm(a); 
end 
end 



if cset==[] 
del_Z=0; 
else 

del_Z=del_K+i*omega(w)*del_C-omega(w)'^2*del_M; 

end 

H_star=Hee-Hec*del_Z*inv(eye(size(Hcc,l)) + Hcc*del_Z)*Hce; 

% Capturing the FRF interested in at each frequency swept thru 
% 

hstar(w)=H_star(frf_index( 1 ),frf_index(2)); 



% Classical method of calculating the FRF by reassembling [Z] and taking the 
% inverse. 

% 

Z_c=K_c+i*omega(w)*C_c-omega(w)'^2*M_c; 

H_c=inv(Z_c); 

H_c=H_c(eset,eset); 
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% Capturing the FRF interested in at each frequency swept thru 
% 

h_c(w)=H_c(frf_index( 1 ),frf_index(2)); 
end 



% Plotting of FRF using Classical & Frequency Synthesis Methods 
% 

setnilogy(ffeq,abs(h_c),'r— ’,ffeq,abs(hstar),’b’) 

title(['Frequency Response Function (H’,int2str(resp_dof),int2str(inp_dof),’)’ ]) 
ylabel(TRF Amplitude (in/lb)') 
xlabelfPrequency (Hz)’) 
legend('Classical',’Synthesis’) 

change_model=inputCWould you like to change your model again? l=yes 2=no '); 



fsynbase.m 



% This program is designed to calculate FRFs by inverting the impedance matrix. 

% Changing the mass, stiffness, and damping matrices, new FRFs will be calculated by 
% reformulating the matrices, and using Ae synthesis method in the frequency domain. 
% 

clear 

% loading of the baseline structure 
% 

dispCSelect the file which contains your [M], [K], and [C] matrices' ) 
pause(2) 

[M_K_C,p]=uigetfile(’*.mat’,’Load [M], [K], [C]’); 

load (M_K_C) 

dofs=l:ndof; 

change_model=l; 
while change_model==l; 

changek=2; changec=2; changem=2; % initializes logic variables to 'no' status 

flag_k=0; flag_c=0; flag_m=0; % initializes change flag 

chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 

cset=0; % initializes cset matrix to zero to avoid null index error in 

% handling the spring to ground if no changes are made to 
% the structure. 



Hey_star=Q; HV_c=[]; % (re)initializes these matrices to an empty mauix 

% to avoid matrix size differences each time the 
% model is changed. 

K_c=K; M_c=M; C_c=C; % initializes change matrices which will be 

% used in the classical method, to the original 
% matrices. 



Kb=zeros(size(K, 1)); % initializes the spring, and damper to base matrices 

Cb=zeros(size(C,l)); % to zero. 



% Designation where base excitation is located, and the spring constant which 
% connects the substructure to the base which is moving. 

% 

bset=input('At what dof(s) are your structure excited? '); 
for b_spring=l:length(bset); 

kb(b_spring)=input(fprintf ('input the value of the spring constant which connects dof 

%g to the base ',bset(b_spring))); 
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end 

cb=zeros( 1 , length (bset)); 



% initializes the damper constant which connects the 
% substructure to the base which is moving, to zero. 



% Correcting the classical [K] & [Kb] matrices to include the spring element 
% connecting the substructure to the base(s). 

% 

for b_spring=l:lengA(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+... 

kb(b_spring); 

Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring)) +... 

kb(b_spring); 

end 

changek=inputCWould you like to change your stiffness matrix? l=yes 2=no '); 
while changek=l 

chk=chk+ 1 ; % counter for determining # change of K matrix 

flag_k=l; 

lk(chk)=input(’input value of added spring (Ibs/in) ’); 

kdofl(chk)=input(’input constrained dof where 1st end of added spring is applied '); 
kdof2(chk)=input('input dof where 2nd end of added spring is applied, b for base, or 

0 for ground '); 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handl^ when asked where the structure is excited. 

% 

if kdofl(chk)~=bset & kdof2(chk)— 'b' 
while kdofl(chk)~=bset & kdof2(chk)=’b’ 

dispCConnecting this dof to the base is changing the basic model.') 
disp(lf this is what you desire, go back and include this dof as an excitation 

point.') 

model_change=inputCWas this change correct? l=yes 2=no ’); 
if model_change==l 

error('Go back and change your basic model.’) 
else 

kdof l(chk)=input(’input constrained dof where 1st end of added damper is 

applied '); 

kdof2(chk)=input(’input dof where 2nd end of added damper is applied, b 

for base, or 0 for ground '); 
end 

end 

end 

% updates the base spring constant and does not count these changes 
% in the cset matrix. 
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% 

if (find(kdofl(chk)=bset))~=0 & kdof2(chk)=='b' 
base=find(kdof 1 (chk)=bset); 
kb(base) = kb(base) + Ik(chk); 
else 

if kdof l(chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdofl(chk)]; % matrix and formation of new c-set matrix, 

end 

if kdof2(chk)~=cset 
cset=[cset kdofZ(chk)]; 
end 
end 

% Classical method of reformulatting the [K] matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk)) + Ik(chk); 
if kdof2(chk)~=0 & kdof2(chk)~=’b’ 
K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)Jcdof2(chk)) - Ik(chk); 
K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk)) - Ik(chk); 
K_c(kdof2(chk)Jcdof2(chk))=K_c(kdof2(chk)Jcdof2(chk)) + Ik(chk); 
end 

% Classical method of reformulatting [Kb] matrix. 

% 

ifkdof2(chk)='b' 

kbb=lk(chk); 

Kb(kdofl(chk),kdofl(chk))=Kb(kdofl(chk),kdofl(chk)) + kbb; 
end 

changek=input(’Are there any other changes to the stiffness matrix? l=yes 2=no '); 
end 

changec=inputCWould you like to change your damping matrix? l=yes 2=no '); 
while changec=l 

chc=chc+ 1 ; % counter for determining # changes to damping matrix. 

flag_c=l; 

lc(chc)=input('input value of added damper (Ibs-secrin) '); 

cdof l(chc)=input(’input constrained dof where 1st end of added damper is applied '); 
cdof2(chc)=input('input dof where 2nd end of added danq)er is applied, b for base, or 
0 for ground '); 

% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handl^ when asked where the structure is excited. 

% 

if cdofl(chc)~=bset & cdof2(chc)=='b' 
while cdofl(chc)~=bset & cdoG(chc)=='b' 

dispCConnecting this dof to the base is changing the basic model.') 
dispCIf this is what you desire, go back and include this dof as an excitation 

point.') 
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model_change=input('Was this change correct? l=yes 2=no '); 
if model_change=l 

error(’Go back and change your basic model.') 

else 

cdofl(chc)=input(’input constrained dof where 1st end of added damper is 

applied ’); 

cdo£2(chc)=input('input dof where 2nd end of added d^per is applied, b 

for base, or 0 for ground '); 
end 
end 
end 



% Updates the base damper constant and does not count these changes 
% in the cset matrix. 

% 

if (find(cdofl(chc)=bset))~=[] & cdof2(chc)='b' 
base=find(cdofl(chc)==bset); 
cb(base) = cb(base) + lc(chc); 
else 

if cdofl(chc)~=cset % comparison of change dof with the c-set 

cset=[cset cdof 1 (chc)]; % matrix and formation of new c-set matrix, 

end 

if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 

% Classical method of reformulatting the [C] matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc)) -i- lc(chc); 
if cdof2(chc)~=0 & cdof2(chc)~='b' 

C_c(cdofl(chc),cdof2(chc))=^_c(cdof l(chc),cdof2(chc)) - lc(chc); 

C_c(cdof2(chc), cdof l(chc))=C_c(cdof2(chc),cdofl (chc)) - lc(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 

% Classical method of reformulatting [Cb] matrix. 

% 

ifcdof2(chc)=’b' 

cbb=lc(chc); 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc)) ebb; 
end 

changec=input(’Are there any other changes to the damping matrix? l=yes 2=no '); 
end 

changem=inputCWould you like to change your mass matrix? l=yes 2=no ’); 
while changem=l 
chm=chm+l; 
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flag_m=l; 

lm(chm)=input('input value of added mass (lbs-sec^2/in) '); 
mdof(chm)=dnput(’input constrained dof where mass is added '); 

% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 

% Classical method of refonnulatting the [M] matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + Im(chm); 

changem=input(’Are there any other changes to the mass matrix? l=yes 2=no ’); 
end 

% Handling of spring added to ground 

% 

ground = find(cset=0); 
cset(ground) = Oi 

% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 

% 

iset=input(’What dofs other than those where change occured, are you interested in? '); 

% Formation of zset (c U b) and eset (i U c U b) matrices. 

% 

zset=[cset bset]; eset=[iset cset bset]; 

% Choosing the FRF interested in 

% 

resp_dof=input('What response dof are you interested in? '); 
ftf_index=find(eset=resp_dof) ; 

% Choosing the type of input and out FRF interested in 

% 

inp_type=inputCIs the base excitation in the form of: l=displacement or 2=acceleration 

resp_type=inputCWhat type of response are you interested in: l=displacement or 

2=acceleration 



% Structural Synthesis method of calculating FRF 
% 



% Calculation of unchanged FRF 
% 

ffeq=.01:l:250+.01; 
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omega=2*pi*freq; 
j=length(omega); 
for w=l:j 

Z=K+i*omega(w)*C-omega(wy^2*M; 

H=inv(Z); 

% Rearrangement of original [H] to form [Heb], [Hez], [Hzz], [Hzb] 

% 

Heb=H(eset,bset); 

Hez=H(eset,zset); 

Hzz=H(zset,zset); 

Hzb=H(zset,bset); 

% Formation of [del_Kc], [del_Cc], [del_Mc], [del_Zc], [del_Zb] 

% & [del_Z] 

% 

del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 

if flag_k=l 
for a=l:chk 

if (find(kdofl(a)==bset))~=D & kdof2(a)=’b’ 

indexk 1 =Q ; % provides for base changes to not be 

indexk2=Q; % included in del_Kc matrix 

else 

indexkl=find(kdofl(a)==cset); % finds where in cset matrix 
indexk2=find(kdof2(a)=cset); % the change in k occurs 
end 

del_Kc(indexkl, indexk l)=del_Kc(indexkl, indexk 1) + lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2) - lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 
end 
end 



if flag_c=l 
for a=l:chc 

if (find(cdofl(a)=bset))~=Q & cdof2(a)=’b’ 
indexcl=[]; 
indexc2=[]; 
else 

indexcl=find(cdofl(a)=cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); % the change in C occurs 
end 

del_Cc(indexcl,indexcl)=del_Cc(indexcl,indexcl) + lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2) - lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2,indexcl) - lc(a); 
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del_Cc(indexc2,indexc2)=del_Cc(indexc2,indexc2) + lc(a); 
end 
end 

if flag_m=l 
fora=l:chm 

indexm=find(mdof(a)=cset); % finds where in cset matrix 

% the change in k occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 
end 
end 

del_Zb=diag(kb+i*omega(w)*cb); % diag is used in the event >1 base 

% excitation is present 

del_Zc=del_Kc+i*omega(w)*del_Cc-omega(w)^2*del_Mc; 

del_Z=[ del_Zc zeros(size(del_Zc,l),size(del_Zb,2)) 
zeros(size(del_Zb, 1 ),size(del_Zc,2)) del_Zb] ; 



% Base excitation Frequency synthesis equation 
% 

Hey_star(:,w)=(Heb*del_Zb -Hez*del_Z*inv(eye(size(Hzz,l))+Hzz*del_Z)*... 
Hzb*del_Zb)*ones(l,length(bset))’; 

if inp_type==l & resp_type=2 

Hey_star(:,w) = Hey_star(:,w)*(-omega(w)'^2); 
elseif inp_type=2 & resp_type— 1 

Hey_star(:,w) = Hey_star(;,w)*(-l/omega(w)'^2); 
end 

% Classical method of calculating the FRF by reassembling [Z] and taking the 
% inverse. 

% 

n=eye(size(M_c, 1 )); 

YV=ones( 1 ,ndof)'; 

Z_c=(K_c+i*omega(w)*C_c-omega(w)'^2*M_c); 

H_c=inv(Z_c); 

HV_c(:,w) = (H_c*(Kb+i*omega(w)*Cb))*YV; 

if inp_type=l & resp_type==2 

HV_c(:,w) = HV_c(:,w)*(-omega(w)^2); 
elseif inp_type==2 & resp_type=l 

HV_c(:,w) = HV_c(:,w)*(-l/omega(w)^2); 
end 
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% Plotting of FRF using Classical & Frequency Synthesis Methods 

% 



semilogy(freq,abs(HV_c(resp_dof,:)),'r"’,freq,abs(Hey_star(frf_index,:)),’b') 
titleCCTrequency Response Function (H',int2str(resp_dof),')' ]) 
if inp_type=l & resp_type=2 
ylabelCFRF Amplitude (g/in)') 
elseif inp_type=2 & resp_type==l 
ylabel(TRF Amplitude (in/g)*) 
else 

ylabelCFRF Amplitude (unitless)') 
end 

xlabelCFrequency (Hz)*) 

legendCClassical', ’Synthesis') 

change_model=inputCWould you like to change your model again? l=yes 2=no ') 
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fsynmain.m 



% This program is designed to calculate the new FRFs of a system after changes in the 
% mass, stiffness, and damping matrices have been made. The new FRFs will be 
% calculated using the syn Aesis method in the frequency domain. 

% 

clear 

% loads the baseline stmcture to be modified, allows for modifications to the baseline 
% structure, and calls the function delta.m 
% 

load change 

% Designation of the frequency range and step size 
% 

% initial step size end 

freq_input = [ .01 1 250+.01 ]; 

freq=freq_input(l):freq_input(2):freq_input(3); 

omega=2*pi*freq; 



[wn,phi,zeta,Mmodal,phi_norm,Cmodal] = modal(M,K,C); 

[hee] = frfmodal(wn,phi,zeta,Mmodal,omega,eset); 

[h_star] = frfsynth(hee,kb,cb,omega,excite,eset,iset,bset,del_Kc,del_Cc,del_Mc); 

% Choosing the type of input and output FRF interested in 
% 

if excite==2 

inp_type=input(ls the base excitation in the form of: l=displacement or 2=acceleration 

resp_type=dnputCWhat type of response are you interested in: l=displacement or 

2=acceleration '); 

ifinp_type=l «& resp_type=2 
h_star = h_star * (-omega(w)'^2); 
elseif inp_type=2 & resp_type=l 
h_star = h_star * (-l/omega(w)'^2); 
end 
end 

resp_dof=inputCWhat response dof are you interested in? '); 

if excite=2 
inp_dof=n; 
else 

inp_dof=inputCWhat input dof are you interested in? ’); 
end 
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[H_star_desired] = frf_sift(eset^esp_dof,inp_dof,excite,h_star); 



semilogy(freq’,abs(H_star_desired),’g’) 

title(['Synthesized Frequency Response Function H',int2str([resp_dof inp_dof]) ]) 
if inp_type=l & resp_type=2 
ylabel(TRF Amplitude (g/in)’) 
elseif inp_type==2 & resp_type=l 
ylabelCF^ Amplitude (in/g)’) 
else 

ylabel(TRF Amplitude (unitless)') 
end 

xlabel(Trequency (rad/s)') 
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% This program provides for the loading of the baseline structure, accepts the user changes 
% to the [Nfl, [K], and [C] matrices, calls the function delta.m, and stores this information 
% to be loaded from another program. 



% 



clear 

changek=2; changec=2; changem=2; % initializes logic variables 
status 



% to 'no' 



flag_k=0; flag_c=0; flag_m=0; 

chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 



cset=0; bset=G; % initializes cset matrix to zero to avoid 

% null index error if no changes are made to 
% the structure and bset to empty matrix to 
% prevent error in the event tlus is not a 
% base excitation and bset does not exist 



kdofl=0; 

cdofl=Q; 

mdof=Q; 

lk=0; Ic^; 

kb0=0; 



kdof2=0; % 

cdof2=D; % 



lm=0, % 

cb0=0; % 

% 



initalizes the matrices which keep 
track of the dofs where changes occur 



initializes the value of added prarmeters 

initializes kb and cb to zero to prevent error in the event 
this is not a base excitation and kb & cb does not exist 



b=num2str('b'); % allows b to be entered as a numerical value, but 

% declares b a string variable to used for 
% comparisons in logic statements. 

dispCAre there already reduced FRF & IRF matrices in the correct format available,') 
dispC'or does the presynthesized FRF &IRF matrices need to be generated using the dof ') 
dispC'locations where you desire to change the structure? ') 

FRF_IRF=input('l=reduced FRF/IRF already exists 2=generate reduced FRF/IRF '); 

% Protects against user making an error in choosing how FRF & IRF is obtained. 

% 

ifFRF_IRF~=[l 2] 
while FRF_IRF~=[1 2 ] 

dispCError in choosing how FRF/IRF is obtained. Choose 1 or 2.') 
FRF_IRF=input('l=reduced FRF/IRF already exists 2=generate reduced FRF/IRF '); 
end 
end 



124 



% Loading the prcsynthesized FRF/IRFs, corresponing frequency vector, and vector of 
% dofs in the eset if reduced FRF/IRF already exists. Loading of M, K, & C matrices if 
% FRF/IRF needs to be generated 
% 

if FRF_IRF=1 

dispCSelect the file which contains your presynthesized FRF/IRFs, a corresponding 
frequency ') 

dispCvector, and the vector of all the dofs that will be in your eset vector.') 
pause(2) 

[hee_dofs_omega,p]=uigetfile('*.mat',’Load [FRF], [IRF], (dofs), {omega}'); 
load (hee_dofs_omega) 

else 

dispCSelect the file which contains your [M], [K], and [C] matrices' ) 
pause(2) 

[M_K_C,p]=uigetfile('*.mat','Load [M], [K], [C]'); 
load (M_K_C) 
dofs=l:ndof; 
end 

excite=inputCHow is the structure excited? l=system 2=base '); 

% Protects against user making an error in choosing the type of excitation. 

% 

ifexcite~=[l 2] 
while excite~=[l 2] 

dispCError in choosing type of excitation. Choose 1 or 2.') 
excite=inputCHow is the stmctuie excited? l=system 2=base '); 
end 
end 

if excite=2 

% Designation where base excitation is located. 

% 

bset=input('At what dof(s) are your structure excited? '); 

% Protects against user making an error in choosing the location of the base excitation. 
% 

for b_dof=l:length(bset) 
if bset(b_dof)~=dofs 

while bset(b_dof)~=dofs 

fprintfCError in choosing location of base excitation at dof %g.\n',... 

bset(b_dof)) 

dispCEither the dof chosen does not exist on your structure, or is not ') 
dispC'included in your list of dofs for the eset vector. ') 
bset(b_dof)=input(fprintf('Choose again the dof for base exitation %g ',... 

b_dof)); 

end 
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end 

end 

% Designation of the spring constant which connects the substructure to the base and 
% protection against user making an error in choosing the value of the constant(s). 

% 

for b_spring=l :length(bset) 

kbO(b_spring)=input(fprintf('input the value of the spring constant which connects dof 

%g to the base ’,bset(b_spring))); 

kb0_zero=find(kb0==0) ; 
if length(kbO)~=b_spring I kbO_zero~=D 

while length(kbO)~=b_spring I kbO_zero~=[] 

fprintf(’You made an error in entering the value for dof %g base 
excitation spring constant.\n',bset(b_spring)) 

kb0(b_spring)=0; 

kbO(b_spring)=input(fprintf(’Re-input the value of the spring constant 
which connects dof %g to the base ’,bset(b_spring))); 
kb0_zero=find(kb0==0); 
end 
end 
end 



cbO=zeros( 1 ,length(bset)); 
end 



% initializes the damper constant which coimects the 
% substructure to the base which is moving, to zero. 



changek=inputCWould you like to make stiffness modifications to the structure? l=yes 

2=no '); 



while changek=l 



% Designation of the spring change and protection against user making an error in 
% choosing the value of the constant(s). 

% 

chk=chk+l; % counter for determining # changes to stiffness matrix. 
lk(chk)=input('input value of added spring (Ib^in) *); 
if length(lk)~=chk 
while length(lk)~=chk 

disp(Tou made an error entering the value for the stiffness change.') 
lk(chk)=0; 

lk(chk)=input('Re-input value of added spring (Ibs/in) ’); 
end 
end 



% User selection of where stiffiiess changes occur. Protects against user making an 
% error in choosing the locations of spring changes. 

% 

kdofl(chk)=input('input dof where 1st end of added spring is applied ’); 
if kdofl(chk)~=dofs 
while kdofl(chk)~=dofs 

dispCError in choosing location of 1st end of added spring. Either the dof 

chosen does') 
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end 

end 



dispCnot exist on your structure, or is not included in your list of dofs for the’) 
dispC'eset vector. ' ) 

kdofl(chk)=input('Re-input constrained dof where 1st end of added spring is 

applied ’); 



kdof2(chk)=input(’input dof where 2nd end of added spring is applied, b for base, or 0 

for ground ’); 

if kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~=’b’ 
while kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~=’b' 

dispCEiror in choosing location of 2nd end of added spring. Either the dof 

chosen does') 

dispCnot exist on your structure, or is not included in your Ust of dofs for the’) 
dispC’eset vector. ’ ) 

kdof2(chk)=input(’Re-input dof where 2nd end of added spring is applied, b for 

base, or 0 for ground ’); 



end 

end 



ifexcite=l & kdof2(chk)=’b’ 
while excite==l & kdof2(chk)=’b’ 

dispCError in choosing location of 2nd end of added spring. You did not 

designate tihis’) 

dispC’as a base excitation structure. Therefore b is not a dof option.’) 
kdoO(chk)=input(’Re-input dof where 2nd end of added spring is applied, b for 

base, or 0 for ground ’); 
end 
end 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset A change from dof to base constitutes an addition of 
% excitation points. This should be handled when asked where the structure is excited. 
% 

if excite==2 & kdofl(chk)~=bset & kdof2(chk)~’b’ 
while kdofl(chk)~=bset & kdof2(chk)=’b’ 

dispC’Connecting this dof to the base is changing the basic model.’) 
dispC’If this is what you desire, go back and include this dof as an excitation 

point.’) 

model_change=inputCWas this change correct? l=yes 2=no ’); 
if model_change==l 

errorC’Go back and change your basic model.’) 
else 

kdofl(chk)=input(’Re-input constrained dof where 1st end of added spring is 

applied ’); 

kdof2(chk)=input(’Re-input dof where 2nd end of added spring is applied, b 

for base, or 0 for ground ’); 
end 
end 
end 
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% Provision for base spring constant changes to not be included in the cset matrix. 

% 

if excite=2 & (find(kdofl(chk)==bset))~=[] & kdof2(chk)==’b’ 
cset=cset; 
else 

if kdof 1 (chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdofl(chk)]; % matrix and formation of new c-set matrix, 

end 

if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 

changek=input('Are there any other stiffiiess modifications to the structure? l=yes 2=no'); 
end 

changec=inputCWould you like to make damping modifications to the structure? l=yes 
2=no '); 

while changec==l 

% Designation of the damper change and protection against user making an error in 
% choosing the value of the constant(s). 

% 

chc=chc+l; % counter for determining # changes to damping matrix. 

lc(chc)=input('input value of added damper (Ibs-secAn) ’); 
if length (lc)~=chc 
while length(lc)~=chc 

dispC'You made an error entering the value for the damper change.') 
lc(chc)=0; 

lc(chc)=input(’Re-input value of added damper Obs-secAn) '); 
end 
end 



% User selection of where damper changes occur. Protects against user making an error 
% in choosing the locations of damper changes. 

% 



cdofl(chc)=input('input constrained dof where 1st end of added damper is applied ’); 
if cdoh(chc)~=dofs 
while cdofl(chc)~=dofs 

dispCError in choosing location of 1st end of added damper. Either the dof 

chosen does') 

dispC'not exist on your structure, or is not included in your list of dofs for the') 
disp('eset vector. ' ) 

cdofl (chc)=input('Re-input constrained dof where 1st end of added damper is 

applied '); 



end 

end 



cdof2(chc)=input('input dof where 2nd end of added damper is applied, b for base, or 0 

for ground '); 
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if cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~=’b’ 
while cdoO(chc^=dofs & cdon(chc)~=0 & cdoK(chc)~='b’ 

dispCError in choosing location of 2nd end of added damper. Either the dof 

chosen does’) 

dispCnot exist on your structure, or is not included in your list of dofs for the') 
dispC’eset vector. ' ) 

cdof2(chc)=input(’Re-input dof where 2nd end of added damper is applied, b for 

base, or 0 for ground ’); 



end 

end 



if excite— 1 & cdof2(chc)='b’ 
while excite==l & cdoO(chc)=='b’ 

dispCError in choosing location of 2nd end of added damper. You did not 

designate this’) 

dispCas a base excitation structure. Therefore b is not a dof option.’) 
cdof2(chc)=inputCRe-input dof where 2nd end of added damper is applied, b for 

base, or 0 for ground ’); 
end 
end 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from these dof to base constitutes an addition 
% of excitation points. This should be handled when asked where the structure is 
% excited. 

% 

if excite=2 & cdofl(chc)~=bset & cdof2(chc)=’b’ 
while cdofl(chc)~=bset & cdof2(chc)==’b’ 

dispCConnecting this dof to the base is changing the basic model.’) 
dispCIf this is what you desire, go back and include this dof as an excitation 

point.’) 

model_change=inputCWas this change correct? l=yes 2=no ’); 
if model_change=l 

error(’Go back and change your basic model.’) 
else 

cdofl(chc)=input(’Re- input constrained dof where 1st end of added damper is 

applied ’); 

cdof2(chc)=input(’Re-input dof where 2nd end of added damper is applied, b 

for base, or 0 for ground ’); 
end 
end 
end 

% Provision for base damper constant changes to not be included in the cset matrix. 

% 

if excite=2 & (find(cdofl(chc)==bset))~=Q & cdof2(chc)=’b’ 
cset=cset; 
else 

if cdof 1 (chc)~=cset % comparison of change dof with the c-set 

cset=[cset cdofl(chc)]; % matrix and formation of new c-set matrix. 
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end 

if cdof2(chc)~=cset 

cset=[cset cdof2(chc)]; 
end 
end 

changec=input('Are there any other damping modifications to the stmcture? l=yes 2=no’); 
end 

changem=inputCWould you like to make a mass modification to the structure? l=yes 
2=no’); 

while changem=l 

% Designation of the mass change and protection against user making an error in 
choosing die 

% value of the constant(s). 

% 

chm=chm+ 1 ; % counter for determining # changes to mass matrix. 

lm(chm)=input(’input value of added mass 0bs-sec^2/in) ’); 
if length(lm)~=chm 
while len^(lm)~=chm 

disp(’You made an error entering the value for the mass change.’) 
lm(chm)=0; 

lm(chm)=inputCRe-input value of added mass (lbs-sec^2/in) ’); 
end 
end 



% User selection of where mass changes occur. Protects against user making an error in 
% choosing the locations of mass changes. 

% 

mdof(chm)=input(’input constrained dof where mass is added '); 
if mdof(chm)~=dofs 
while mdof(chm)~=dofs 

dispCError in choosing location of mass change. Either the dof chosen does 

not’) 

disp(’exist on your stracture, or is not included in your list of dofs for the’) 
dispCeset vector. ’ ) 

mdof(chm)=input(’input constrained dof where mass is added ’); 
end 
end 

if mdof(chm)~=cset % comparison of change dof with the c-set matrix 

cset=[cset mdof(chm)]; % and formation of new c-set matrix, 

end 

changem=input(’Are there any other mass modifications to the stmcture? l=yes 2=no ’); 
end 

% Handling of spring added to ground to eliminate zero from cset 
% 
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ground = find(cset==0); 
cset(ground) = []; 

dispCThe following is a list of dofs were you have made changes to your structure.’) 

disp(This represents your {cset} vector:') 

disp(cset) 

% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 

% 

ifFRF_IRF=l 

dispCSince you inputed your own FRF/IRF matrices, your iset is chosen as all dofs 
remaining ') 

dispCin your ese^dof vector after extracting the cset vector.') 
iset=dofs; 
for a=l :length(cset) 
extract(a)=find(cset(a)=dofs); 
end 

iset(extiact)=D; 

else 

iset=input( What dofs other than those where change occured, are you interested in? '); 
end 

% Formation of zset (c U b) and eset (i U c U b) matrices. 

% 

zset=[cset bset]; eset=[iset cset bset]; 

[del_Kc,del_Cc,del_McJcb,cb] = 

delta(cset,bset,kdof 1 4cdof2,lk,cdof 1 ,cdof2,lc,mdof,lm,kb0,cb0); 

save change M K C ndof del_Kc del_Cc del_Mc kb cb iset cset bset zset eset excite 
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modal.m 



function [wn,phi,zeta>lmodal,phi_nonn,Cmodal] = modal(M,K,C) 

% This function is designed to calculate the natural frequencies and 
% mode shapes 
% 

ndof=size(M,l); 

[PHI,lam]=eig(M\K); 
for j=l;ndof; 

lambda(j)=lam(jj); 

end 

[wn,I]=sort(sqrt(lambda)); 
wn=rot90(wn,- 1 ); 
for j=l;ndof; 

phi(:,j)=PHI(:,ia)); 

end 

% Formulation of modal mass and damping matrices 
% 

Mmodal=phi'*M*phi; 
phi_norm=phi/(Mmodal'^. 5); 

Mmodal=diag(l^odal); 

Cmodal=phi'*C*phi; 

Cmodal=^ag(Cmodtd); 

%zeta0=0.0; 

%zeta=zetaO*ones( 1 ,length(Cmodal)); 
zeta=Cmodal./(2*Mmodal.*wn); 



% Constant modal damping ratio 

% Modal damping ratio derived 
% from proportional damping matrix [C] 
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frfmodal.m 



function [h] = fffmodal(wn,phi,zeta,Mmodal,omega,frf_dof) 

% This function is designed to calculate the FRF using mode shapes and 
% summing the FRF over a number of modes. 

% 



ncol=(length(fTf_dof)*(length(fTf_dof)+l))/2; 

% Calculation of frequency response function 
% 

h=zerosGength(omega),ncol); 
triangle matrix. 



% Initialization to zero, and sizing 

% of the freq x lower 



h_old=0; h_new=0; % initialization of the matrices to 

% be used to determine when modal 
% convergence has occured. 

h_check=l; % Initialization of the stop criteria for summing the 

% modes 

nmodes=size(phi); b=0; % Defining the number of modes that exist and 

% initializing b which keeps track of the 
% of modes used for convergence of the FRF. 

HOdprime_aa=zeros(length(fif_dof)); 

for b=l:nmodes 

% Elimination of unnessecary elements and rearrangement of original 
% mode shapes {phi) to form {phi_red} and [num_red] which is now an 
% frf_dof X firf_dof matrix. 

% 

num_red=phi(firf_dof,b)*phi(frf_dof,b)'; 
for w=l :length(omega) 

den=(wn(by^2 - omega(wy^2 + 2*i*zeta(b)*wn(b)*omega(w))*Mmodal(b); 
H_red=num_red./den; 

% Changes the [H_red] FRFs from an frf_dof x frf_dof matrix to 
% a freq x lower triangle matrix. 

% 

H_lowtri=tril(H_red); % pulls out the lower triangle 

% and places zeros everywhere else 

symtry = find(H_lowtri==0); % finds zero values in the 
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% lower triangle matrix 



H_lowtri(symtry) = []; 

h_mode(w,:)=H_lowtri; 

end 

h = h + h_mode; 
end 



% deletes all values of zero and 
% turns the remaining elements 
% into a 1 X length(lower triangle) 
% vector 

% saves the lower triangle vector 
% at each freq. These raFs are for 
% a single mode only. 



% Summing of each modal FRF 



134 



frfsynth.m 



function [h_star] = frfsynth(heeJkb,cb,omega,excite,eset,iset,bset,... 
del_Kc,del_Cc,del_Mc) 

% This fimction is designed to calculate the new FEtF using the synthesis 
% method. 

% 

for w=l :length(omega) 

count=0; % initializes the counter which keeps track of 

% which column of the hee matrix is to be 
% placed into the refomaticMi of matrix [Hee] 

for col = l:length(eset) % reforms the [Hee] lower 

for row = col:length(eset) % triangle matrix 

count=count + 1; 

Hee_lowtri(row,col)=hee(w,count); 

end 

end 

% reforms the original symmetric [Hee] matrix 
% 

Hee=Hee_lowtri; 

[sym_row,sym_col] = find(Hee_lowtri==0); 
for n=l :length(sym_row) 

Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 

end 

e=l:length(eset); 

c=length(iset)+ 1 :length(eset)-length(bset); 
b=length(eset)-length(bset)+l:length(eset); 
z=length(iset)+l :length(eset); 

% Rearrangement of [Hee] to form [Hec], [Hcc], [Hee] if the structure 
% is system excited. 

% 

if excite=l 
Hec=Hee(e,c); 

Hce=Hee(c,e); 

Hcc=Hee(c,c); 

end 

% Rearrangement of [Hee] to form [Heb], [Hez], [Hzz], [Hzb] if the 
% structure is base excited. 

% 

if excite=2 
Heb=Hee(e,b); 
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Hez=Hee(e,z); 

Hzz=Hee(z,z); 

Hzb=Hee(z,b); 

end 



% Formation of [del_Zc], [del_Zb] & [del_Z] 

% 

del_Zc=del_Kc+i*omega(w)*del_Cc-omega(w)'^2*del_Mc; 
if excite==2 

del_Zb=diag(kb+i*omega(w)*cb); % diag is used in the event 

%>1 base excitation 
% is present 

del_Z=[ del_Zc zeros(size(del_Zc,l),size(del_Zb,2)) 
zeros(size(del_Zb,l),size(del_Zc,2)) del_Zb]; 
end 

% Determining which synthesis equations are to be used 
% 

if excite==l 



% Structure excitation Frequency synthesis equation 
% 

H_star=Hee-Hec*del_Zc*inv(eye(size(Hcc,l)) + Hcc*del_Zc)*Hce; 
else 

% Base excitation Frequency synthesis equation 
% 

h_star(:,w)=(Heb*del_Zb - 

Hez*del_Z*inv(eye(size(Hzz, 1 ))+Hzz*del_Z)*Hzb*del_Zb)*ones( 1 ,length(bset))'; 
end 

if excite=l 

% Changes the [H*] FRFs from an eset x eset matrix to a freq x lower 
% triangle matrix, and creates an index matrix which keeps track of 
% the original coordinates of the FRFs in the [H*] matrix. 

% 

Hstar_lowtri=tril(H_star); % pulls out the lower triangle and places 

% zeros everywhere else 

symtry = find(Hstar_lowtri=0); % finds zero values in the lower 

% triangle matrix 

Hstar_lowtri(symtry) = []; % deletes all values of zero and turns the 
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h_star(w,:)=Hstar_lowtri; 

end 

end 



% remaining elements into a 1 x length(lower 
% triangle) vector 

% saves the lower triangle vector at each freq 



if excite=l 

% changes h_star matrix so ffeqs are in each column and the FRFs are 
% in the rows 
% 

h_star=flipud(rot90(h_star)); 

end 
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frf sift.m 



function [H_star_desired] = fif_sift(eset,resp_dof,inp_dof,excite,h_star) 



% Locates the positions in the eset that the response and input dofs are referring. 

% These locations correspond to the matrix location row and column of the desired 
% FRF in the [H*] matrix which is eset x eset and partitioned. Or these locations 
% refer to the row of the [H*] matrix which is eset x 1 for a base excitation. 

% These indices sift through the [h*] matrix rows to find the desired FRF. 

% 



% Creates an index matrix which keeps track of the original coordinates of the FRFs 
% in the [H_star] matrix. 

% 

count=0; % initializes the counter which keeps track of 

% which column of the h matrix the lower triangle [H_red] 

% FRF went into 



for col = 1 :length(eset) % creates a matrix which matches up 

for row = col:length(eset) % each element of the lower triangle 

count=count +1; % H_star matrix with the rows in the 

hstar_index(count,:)=[row col]; % h_star matrix that they were placed in. 

end 
end 

fif_index=find(eset==resp_dof I eset=inp_dof); 

% prevents index error in the event response dof and input dof are the same 
% 

if resp_dof=inp_dof 
fif_index=[fif_index firf_index]; 
end 

if excite=l 

% Uses hstar_index and finds the row in the lower triangle x freq matrix that the 
% desired FRF is located. 

% 

fif_row=find((frf_index(l)=hstar_index(:,l) & ftf_index(2)=hstar_index(:,2)) I 
(frf_index(l)=hstar_index(:,2) & fif_index(2)=hstar_index(:,l))); 

H_star_desired = h_star(fif_row,:); 

else 

H_star_desired = h_star(frf_index,:); 
end 
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APPENDIX B. STATIC DISPLACEMENT SYNTHESIS COMPUTER 



CODES 



The following is a list and a brief description of the main MATLAB 

computer codes that were written in order to perform the static displacement 

s)mthesis calculations. 

• Guyan xstat.m - uses Guyan reduction methods to calculate the static 
displacements on a structure. 

• ssynmain.m - the main program which calls the modularized static 
displacement synthesis programs and does the post processing of the 
results. 

• statchange.m - loads the baseline structiore to be modified, allows for 
modifications to the baseline structure, and calls the function statdelta.m 

• statdelta.m - forms the modification matrices which will be used to form 
impedance matrices. 

• modal.m - solves for the structure's natural frequencies, mode shapes and 
other modal matrices and vectors. 

• frf modal.m - calculates the structure's FRFs. 

• statsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized static displacements. 

The full codes are contained on subsequent pages. Any codes that were 

previously presented will not be repeated. 
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Guyanxstat.m 



% The purpose in this program is to use Guyan reduction methods in order to calculate the 
% static sisplacements on a stmcture. 

% 

clear 

load fe_add 

% Stiffness changes to the stmcture 
% 

K(3,3)=K(3,3)-1000; 

K(93,93)=K(93,93)-1000; 

K(3,3)=K(3,3)+le4; 

K(93,93)=K(93,93)+le4; 

K(66,66)=K(66,66)+500; 

K(66, 109)=K(66, 1 09)-500; 
tic 

g=ones(size(M,l),l)*(-386.4); 
for rot=l:3:size(M,l)-3 
g(rot)=0; 
end 

for rot=2:3:size(M,l)-2 
g(rot)=0; 
end 

F=M*g; 
dofset=l:ndof; 
oset=dofset; 

aset=[3 18 93 108 66 109]; 
for a=l:length(aset) 
ldll(a)=find(aset(a)==dofset) ; 
end 

osetOdll)=[]; 
nset=[aset oset]; 

Knn=K(nset,nset); Fnn=F(nset); 

Kaa=K(aset,aset); Kao=K(aset,oset); Koo=K(oset,oset); Koa=K(oset,aset); 
Toa=(- 1 ) *inv(Koo)*Koa; 

T=[eye(length(aset));Toa]; 

Kred_Guyan=T'*Knn*T; Fred_Guyan=T’*Fnn; 

xstat_Guyan=Kred_Guyan\Fred_Guyan; 



K(18,18)=K(18,18)-1000; 

K(108,108)=K(108,108)-1000; 

K(18,18)=K(18,18)+le4; 

K(108,108)=K(108,108)+le4; 

K(109,109)=K(109,109)+500; 

K(109,66)=K(109,66)-500; 
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ssynmain.m 



% This program is designed to calculate the new static displacements of a system after 
% changes in the mass, stiffness, and damping matrices have been made. Tlie new static 
% displacements will be calculated using the static displacement synthesis method 
% 

clear 

load statchange 
stat_omega=. 1 : . 1 : .4; 

C=zeros(size(C)); 

[stat_wn,stat_phi,stat_zeta,stat_Mmodal] = modal(M,K,C); 
tO = clock; 

[stat_hee] = fifmodal(stat_wn,stat_phi,stat_zeta,stat_Mmodal,stat_omega,eset); 



[z_stat,Fred,Kred,Mred,HstarO_dp] = statsynth(stat_hee,kb,stat_omega,eset,iset,bset,... 

del_Kc,del_Mc); 

tl = clock; 

stat_syn_time=etime(tl ,t0) 
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statchange.m 



% This program provides for the loading of the baseline structure, accepts the user changes 
% to the [Nf|, [K], and [C] matrices, calls the function statdelta.m, and stores this 
% information to be loaded from another program. 



% 



clear 

changek=2; changec=2; changem=2; % initializes logic variables 
status 



% to ’no’ 



flag_k=0; flag_c=0; flag_m=0; 

chk=0; chc=0; chm=0; % initializes counter for # of changes 

% to matrices 



cset=0; 


bset=D; 


% initializes cset matrix to zero to avoid 
% null index error if no changes are made to 
% the stmcture and bset to empty matrix to 
% prevent error in the event this is not a 
% base excitation and bset does not exist 


kdofl=D; 

cdofl=D; 

mdof=[]; 


kdof2=D; 

cdof2=Q; 


% initalizes the matrices which keep 
% track of the dofs where changes occur 


lk=0; lc=0; 


lm=0; 


% initializes the value of added praimeters 


kb0=0; 


cb=0; 


% initializes kb and cb to zero to 
% prevent error in the event this is not a 
% base excitation and kb & cb does not exist 



b=num2str(’b’); % allows b to be entered as a numerical value, but 

% declares b a string variable to used for 
% comparisons in logic statements. 

% Used to set the range of possible dofs for the user to choose when making changes 
% to his structure. 

% 

dispC’Is there already a reduced FRF matrix in the correct format available,’) 
disp(’or does the presynthesized FRF matrix need to be generated using the dof ’) 
disp(’locations where you desire to change the structure? ’) 

FRF=input(’l=reduced FRF already exists 2=generate reduced FRF ’); 

% Protects against user making an error in choosing how FRF is obtained. 
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% 

if FRF~=[1 2] 
while FRF~=[1 2] 

dispCError in choosing how FRF is obtained. Choose 1 or 2.') 
FRF=input('l=reduced FRF already exists 2=generate reduced FRF ’); 

end 
end 

% Loading the presynthesized FRFs, corresponing frequency vector, and vector of dofs in 
% the eset if reduced FRF already exists. Loading of M, K, & C matrices if FRF needs to 
% generated 
% 

ifFRF=l 

dispC’Select the file which contains your presynthesized FRFs, a corresponding 
frequency ’) 

disp(Vector, and the vector of all the dofs that will be in your eset vector.’) 
pause(2) 

[hee_dofs_omega,p]=uigetfile(’*. mat’, 'Load [FRF], {dofs}, {omega}’); 
load (hee_dofs_omega) 

else 

dispC’Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 

[M_K_C,p]=uigetfile(’*.mat’,’Load [M], [K], [C]’); 
load (M_K_C) 
dofs=l:ndof; 
end 

excite=input(’How is the stmcture excited? l=system 2=base ’); 

% Protects against user making an error in choosing the type of excitation. 

% 

if excite~=[l 2] 
while excite~=[l 2] 

dispCError in choosing type of excitation. Choose 1 or 2.’) 
excite=inputCHow is Ae structure excited? l=system 2=base ’); 
end 
end 

if excite=2 

% Designation where base excitation is located. 

% 

bset=input(’At what dof(s) are your structure excited? ’); 

% bset=l; 

% Protects against user making an error in choosing the location of the base excitation. 

% 

for b_dof=l:length(bset) 
if bset(b_dof)~=dofs 

while bset(b_doO~=dofs 

fprintfCError in choosing location of base excitation at dof %g.\n’,... 
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bset(b_dof)) 

dispCEither the dof chosen does not exist on your structure, or is not ') 
dispCincluded in your list of dofs for the eset vector. ') 
bset(b_dof)=input(fprintf('Choose again the dof for base exitation %g 

b_dof)); 



end 



end 

end 



% Designation of the spring constant which connects the substructure to the base which 
is moving. 

% and protection against user making an errcr in choosing the value of the constant(s). 

% 

for b_spring=l:length(bset) 

kbO(b_spring)=input(fprintf('input the value of the spring constant which connects dof 

%g to the base ’,bset(b_spring))); 

kb0_zero=find(kb0=0) ; 
if length(kbO)~=b_spring I kbO_zero~=Q 

while length(kbO)~=b_spring I kbO_zero~=Q 

fprintf(’You made an error in entering the value for dof %g base 
excitation spring constant.\n’,bset(b_spring)) 

kb0(b_spring)=0; 

kbO(b_spring)=input(fprintf(’Re-input the value of the spring constant 
which connects dof %g to the base ’,bset(b_spring))); 
kb0_zero=find(kb0==0); 
end 

end 

end 



cb=zeros( 1 ,length(bset)); 
end 



% initializes the damper constant which connects the 
% substmcture to the base which is moving, to zero. 



changek=inputCWould you like to make stiffness modifications to the structure? l=yes 
2=no ’); 

while changek— 1 

% Designation of the spring change and protection against user making an error in 
choosing the 

% value of the constant(s). 

% 

chk=chk+l; % counter for determining # changes to stiffness matrix. 
lk(chk)=dnput(’input value of added spring (Ibs/in) '); 
if length(lk)~=chk 
while length (lk)~=chk 

disp(Tou made an error entering the value for the stiffness change.’) 
lk(chk)=0; 

lk(chk)=input(’Re-input value of added spring (Ibsfin) *); 

end 

end 
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% User selection of where stiffness changes occur. Protects against user making an 
%error in choosing 
% the locations of spring changes. 

% 



kdofl(chk)=input(’input constrained dof where 1st end of added spring is applied ’); 
if kdoh(chk)~=dofs 
while kdofl(chk)~=dofs 

dispCError in choosing location of 1st end of added spring. Either the dof 

chosen does') 

dispCnot exist on your structure, or is not included in your list of dofs for the’) 
disp('eset vector. ' ) 

kdofl(chk)=inputCRe-input constrained dof where 1st end of added spring is 

applied '); 



end 

end 



kdof2(chk)=input('input dof where 2nd end of added spring is applied, b for base, or 0 
for ground '); 

if kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~='b’ 
while kdof2(chk)~=dofs & kdof2(chk)~=0 & kdof2(chk)~='b’ 

dispCEiror in choosing location of 2nd end of added spring. Either the dof 

chosen does’) 

disp(’not exist on your structure, or is not included in your list of dofs for the’) 
dispC’eset vector. ’ ) 

kdof2(chk)=input(’Re-input dof where 2nd end of added spring is applied, b for 

base, or 0 for ground ’); 



end 

end 



ifexcite=l &kdof2(chk)==’b’ 
while excite==l & kdof2(chk)=’b’ 

disp(Error in choosing location of 2nd end of added spring. You did not 

designate this’) 

dispCas a base excitation structure. Therefore b is not a dof option.’) 
kdof2(chk)=input(’Re-input dof where 2nd end of added spring is applied, b for 

base, or 0 for ground ’); 
end 
end 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset A change from dof to base constitutes an addition of 
% excitation points. This should be handled when asked where the structure is excited. 
% 

if excite==2 & kdofl(chk)~=bset & kdof2(chk)=='b’ 
while kdofl(chk)~=bset & kdof2(chk)='b’ 

dispCConnecting this dof to the base is changing the basic model.’) 
dispCIf this is what you desire, go back and include this dof as an excitation 

point.’) 

model_change=inputCWas this change correct? l=yes 2=no ’); 
if model_change=l 
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error('Go back and change your basic model.') 
else 

kdofl(chk)=input(’Re-input constrained dof where 1st end of added spring is 

applied ’); 

kdof2(chk)=input('Re-input dof where 2nd end of added spring is applied, b 

for base, or 0 for ground '); 

end 

end 

end 

% Provision for base spring constant changes to not be included in the cset matrix. 

% 

if excite=2 & (find(kdofl(chk)==bset))~=[] & kdof2(chk)=='b’ 
cset=cset; 
else 

if kdofl(chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdofl(chk)]; % matrix and formation of new c-set matrix, 

end 

if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 
end 

changek=input(’Are there any other stiffness modifications to the structure? l=yes 2=no 
end 

changec=inputCWould you like to make damping modifications to the structure? l=yes 
2=no '); 

while changec=l 

% Designation of the damper change and protection against user making an error in 
choosing tiie 

% value of the constant(s). 

% 

chc=chc+ 1 ; % counter for determining # changes to damping matrix. 

lc(chc)=input('input value of added damper (lbs- sec/in) '); 
if length(lc)~=chc 
while length(lc)~=chc 

dispCYou made an error entering the value for the damper change.') 
lc(chc)=0; 

lc(chc)=input('Re-input value of added damper Obs-sec^) '); 
end 
end 



% User selection of where damper changes occur. Protects against user making an error 
% in choosing the locations of damper changes. 

% 

cdofl(chc)=input('input constrained dof where 1st end of added damper is applied '); 
if cdoh(chc)~=dofs 
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while cdofl(chc)~=dofs 

dispCError in choosing location of 1st end of added damper. Either the dof 
chosen does') 

dispC’not exist on your structure, or is not included in your list of dofs for the') 
dispC'eset vector. ' ) 

cdofl(chc)=input('Re-input constrained dof where 1st end of added damper is 
applied '); 
end 
end 

cdof2(chc)=input('input dof where 2nd end of added damper is applied, b for base, or 0 
for ground '); 

if cdof2(chc)~=dofs & cdof2(chc)~=0 & cdof2(chc)~='b' 
while cdoO(chc)~=dofs & cdoO(chc)~=0 & cdoC(chc)~='b' 

dispCError in choosing location of 2nd end of added damper. Either the dof 
chosen does') 

dispCnot exist on your structure, or is not included in your list of dofs for the') 
dispCeset vector. ' ) 

cdof2(chc)=input('Re-input dof where 2nd end of added damper is applied, b for 
base, or 0 for ground '); 
end 
end 

if excite=l & cdof2(chc)='b' 
while excite==l & cdoO(chc)=='b' 

dispCError in choosing location of 2nd end of added damper. You did not 
designate this') 

dispCas a base excitation structure. Therefore b is not a dof option.') 
cdof2(chc)=input('Re-input dof where 2nd end of added damper is applied, b for 
base, or 0 for ground '); 
end 
end 

% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from these dof to base constitutes an addition 
of 

% excitation points. This should be handled when asked where the stmcture is excited. 

% 

if excite=2 & cdofl(chc)~=bset & cdof2(chc)=='b' 
while cdofl(chc)~=bset & cdof2(chc)=='b' 

dispCConnecting this dof to the base is changing the basic model.') 
dispC'If this is what you desire, go back and include this dof as an excitation 

point.') 

model_change=inputCWas this change correct? l=yes 2=no '); 
if model_change=l 

error('Go back and change your basic model.') 
else 

cdofl(chc)=input('Re-input constrained dof where 1st end of added damper is 

applied '); 
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cdof2(chc)=input(’Re-input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground '); 
end 
end 
end 

% Provision for base damper constant changes to not be included in the cset matrix. 

% 

if excite=2 & (find(cdofl(chc)==bset))~=Q & cdof2(chc)=’b’ 
cset=cset; 
else 

if cdof 1 (chc)~=cset % comparison of change dof with the c-set 

cset=[cset cdofl(chc)]; % matrix and formation of new c-set matrix, 

end 

if cdof2(chc)~=cset 

cset=[cset cdof2(chc)]; 
end 
end 

changec=input('Are there any other damping modifications to the structure? l=yes 2=no 
end 

changem=inputCWould you Wee to make a mass modification to the structure? l=yes 2=no 
while changem=l 

% Designation of the mass change and protection against user making an error in 
choosing die 

% value of the constant(s). 

% 

chm=chm+l ; % counter for determining # changes to mass matrix. 

lm(chm)=input(’input value of added mass (lbs-sec'^2/in) '); 
if length(lm)~=chm 
while len^(lm)~=chm 

dispCYou made an error entering the value for the mass change.') 
lm(chm)=0; 

lm(chm)=input('Re-input value of added mass (lbs-sec'^2/in) ’); 
end 
end 



% User selection of where mass changes occur. Protects against user making an error in 
% choosing the locations of mass changes. 

% 

mdof(chm)=input('input constrained dof where mass is added ’); 
if mdof(chm)~=dofs 
while mdof(chm)~=dofs 

dispCError in choosing location of mass change. Either the dof chosen does 

not') 

dispCexist on your structure, or is not included in your list of dofs for the') 
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dispCeset vector. ’ ) 

mdof(chm)=input('input constrained dof where mass is added '); 



end 



end 



if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 



% comparison of change dof with the c-set matrix 
% and formation of new c-set matrix. 



changem=input('Are there any other mass modifications to the structure? l=yes 2=no ’); 
end 

% Handling of spring added to ground to eliminate zero finom cset 
% 

ground = find(cset=0); 
cset(ground) = []; 

dispCThe following is a list of dofs were you have made changes to your structure.’) 

disp(This represents your {cset} vector:’) 

disp(cset) 

% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 



disp(’Since you inputed your own FRF matrix, your iset is chosen as all dofs remaining 



disp(’in your eset/dof vector after extracting the cset vector.’) 
iset^ofs; 
for a=l:length(cset) 
extract(a)=find(cset(a)=dofs); 
end 

iset(extraa)=Q; 

else 

iset=input(’What dofs other than those where change occured, are you interested in? ’); 
end 



% 

ifFRF=l 



’) 



% Formation of zset (c 
% 

zset=[cset bset]; 



U b) and eset (i U c U b) matrices. 
eset=[iset cset bset]; 



[del_Kc,del_Mc4cb] = statdelta(cset,bset,kdofl,kdof2,lk,mdof,lm,kb0); 
save statchange M K C ndof del_Kc del_Mc kb iset cset bset zset eset excite 
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statdelta.m 



% The purpose of this function is to form the modification matrices 
% 



function [del_Kc,del_McJcb] = statdelta(cset,bset,kdofljcdof2,lk,mdof,lmjcb0) 



% Formation of [del_Kc], [del_Cc], & [del_Mc] and updating of kb & cb 
% 



del_Kc=zeros(length(cset)); 

del_Mc=zerosGength(cset)); 

matrix 



% intializes the change matrices 
% to zero and sets their sizes 

% as [cset X cset] 



kb=kbO; % initializes the base excitation spring constant 

% to the initial value inputed when structure was 
% formed. 



if kdofl~=D 

for a=l:length(kdofl) 

if (find(kdofl(a)=bset))~=n & kdof2(a)=='b* 
base=find(kdof 1 (a)==bset); 
kb(base) = kb(base) + lk(a); 



else 

indexk 1 =find(kdof 1 (a)==cset); % finds where in cset matrix 
indexk2=find^dof2(a)=cset); % the change in k occurs 

del_Kc(indexkl,indexkl)=del_Kc(indexkl,indexkl) + lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2) - lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 

end 

end 

end 



if mdof~=[] 
for a=l:length(mdof) 

indexm=find(mdof(a)==cset); % finds where in cset matrix 

% the change in M occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 
end 
end 
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statsynth.m 



function [z_stat,Fred,Kred,Mred,HstarO_dp] = statsynth(hee,kb,omega,eset,iset,bset,.„ 

del_Kc,del_Mc) 

% This function is designed to calculate the new static displacement using the synthesis 
% method. 

% 



pull=D; 



for w=l:length(omega) 

count=0; % initializes the counter which keeps track of 

% which column of the hee matrix is to be 
% placed into the refomation of matrix [Hee] 

for col = l:length(eset) % reforms the [Hee] lower 

for row = col:length(eset) % triangle matrix 

count=count + 1; 

Hee_lowtri(row,col)=hee(w,count); 

end 

end 

% reforms the original symmetric [Hee] matrix 
% 

Hee=Hee_lowtri; 

[sym_row,sym_col] = find(Hee_lowtri==0); 
for n=l:len^(sym_row) 

Hee(sym_row(n),sym_col(n))=Hee(sym_col(n),sym_row(n)); 

end 



e=l:length(eset); 

c=length(iset)+ 1 :length(eset); %-length(bset) ; 

% b=length(eset)-length(bset)+ 1 :length(eset); 

% z=length(iset)+l:lengA(eset); 

% Rearrangement of [Hee] to form [Hec], [Hcc], [Hee] if the structure 
% is system excited. 

% 

% if excite=l 
Hec=Hee(e,c); 

Hce=Hee(c,e); 

Hcc=Hee(c,c); 

% end 

% Formation of [del_Zc], [del_Zb] & [del_Z] 

% 
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del_Zc=del_Kc-omega(w)^2*del_Mc; 

del_Zb=diag(kb); % diag is used in the event >1 base excitation 

% is present 

del_Z=[ del_Zc zeros(size(del_Zc,l),size(del_Zb,2)) 
zeros(size(del_Zb, l),size(del_Zc,2)) del_^]; 



% Structure excitation Frequency synthesis equation 
% 

H_star=Hee-Hec*del_Z*inv(eye(size(Hcc,l)) + Hcc*del_Z)*Hce; 

% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 
% eset, and these rows and columns are deleted from H_star. This is done 
% because the redundant dof column(s) and row(s) in H_star creates a singular 
% matrix which is not invertable and Kred cannot be determined. 

% 

for u=l .'length (bset) 

locate=find(bset(u)==eset); 
if length(locate) > 1 
pull=[pull locate(2)]; 
end 

end 

H_star(pull,:)=n; 

H_star(:,pull)=D; 

if w=l 

Kred=inv(H_star); 

H0=H_star; 

end 

ifw=2 

Hl=H_star, 

end 

ifw=3 

H2=H_star; 

end 

if w=4 
H3=H_star; 
end 

Hstar0_dp=(-H3+4*H2-5*Hl+2*H0)./(omega(2)-omega(l))''2; 

Mred=.5*Kred*HstarO_dp*Kred; 



ga=ones(size(Mred,l), l)*(-386.4); 

Fred=Mred*ga; 

z_stat=Kred^Fred; 
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APPENDIX C TIME DOMAIN SYNTHESIS COMPUTER CODES 



The following is a list and a brief description of the main MATLAB 

computer codes that were written in order to perform the time domain 

synthesis calculations. 

• tsvnconv.m - performs the time domain synthesis on a structure which 
experiences a base excitation. The dynamic responses are compared to the 
responses calculated using the classical convolution integral method. 

• tsynode45.m - performs the time domain synthesis on a structure which 
experiences a base excitation. The dynamic responses are compared to the 
responses calculated using the classical direct integration method. 

• tsynmain.m - the main program which calls the modularized time 
synthesis programs and does the post processing of the results. 

• change.m - loads the baseline structure to be modified, allows for 
modifications to the baseline structure, and calls the function delta.m 

• delta.m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

• modal.m - solves for the structure's natural frequencies, mode shapes and 
other modal matrices and vectors. 

• impmodal.m - calculates the structure's IRFs. 

• buildA.m - builds the quadrature matrices. 
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• fBlastForcing.m - inputs the base excitation as a function of time. 

• timesynth.m - performs the S3mthesis on the baseline structure and 
returns the synthesized transient responses. 

• time sift.m - selects the dynamic response that the user wants to perform 
post processing on. 

The full codes are contained on subsequent pages. Any codes that were 
previously presented will not be repeated. 
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tsynconv.m 



% The purpose of this program is to perform the time domain synthesis on a structure and 
% compare the dynamic responses to those calculated using a convolution integral method 
% 

clear 

dispCSelect the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 

[M_K_C,p]=uigetfile(’*.mat’,'Load [M], [K], [C]’); 
load (M_K_C) 
dofs=l:ndof; 

C=00*C; 



[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,C); 

changek=2; changec=2; changem=2; % initializes logic variables 

% to 'no' status 



flag_k=0; flag_c=0; flag_m=0; 

chk=0; chc=0; chm=0; % initializes counter for# of changes 

% to matrices 



cset=0; % initializes cset matrix to zero to avoid null index error in 

% handling the spring to ground if no changes are made to 
% the structure. 



Hey_star=D; HV_c=[]; 



K_c=K; M_c=M; C_c=C; 



Kb=zeros(size(K, 1)); 
Cb=zeros(size(C,l)); 



% (re)initializes these matrices to an empty matrix 
% to avoid matrix size differences each time the 
% model is changed. 

% initializes change matrices which will be used in 
% the classical method, to the original matrices. 

% initializes the spring, and damper to base matrices 
% to zero. 



K_vec = zeros(ndof,l); 

b=num2str(’b'); % allows b to be entered as a numerical value, but 

% declares b a string variable to used for 
% comparisons in logic statements. 

clear M K C 



% Designation where base excitation is located, and the spring constant which 
% connects the substructure to the base which is moving. 
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% 

bset=input('At what dof(s) are your structure excited? '); 
for b_spring=l:length(bset); 

kb(b_spring)=input(fprintf('input the value of the spring constant which connects dof 

%g to the base ’,bset(b_spring))); 
end 

cb=zeros(l ,length(bset)); % initializes the damper constant which connects the 

% substmcture to the base which is moving, to zero. 

% Correcting the classical [K] & [Kb] matrices to include the spring element 
% connecting the substructure to the base(s). 

% 

for b_spring=l:length(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 
Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring)) + 

kb(b_spring); 

K_vec(bset(b_spring))=K_vec(bset(b_spring)) + kb(b_spring); 
end 

changek=inputCWould you like to change your stiffness matrix? l=yes 2=no '); 
while changek==l 

chk=chk+ 1 ; % counter for determining # changes 

flag_k=l; 

lk(chk)=input(’input value of added spring (IbsAn) '); 

kdofl(chk)=input(’input constrained dof where 1st end of added spring is applied ’); 
kdof2(chk)=input('input dof where 2nd end of added spring is applied, b for base, or 

0 for ground ’)» 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handl^ when asked where the stmcture is excited. 

% 

if kdofl(chk)~=bset & kdof2(chk)=’b’ 

while kdofl(chk)~=bset & kdof2(chk)==’b’ 
dispCConnecting this dof to the base is changing the basic model.’) 
disp(’If this is what you desire, go back and include this dof as an excitation 

point.’) 

model_change=input(’Was this change correct? l=yes 2=no ’); 
if model_change=l 

error(’Go back and change your basic model.’) 

else 

kdofl(chk)=input(’input constrained dof where 1st end of added 
damper is applied ’); 

kdof2(chk)=input(’input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground ’); 
end 

end 

end 
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% updates the base spring constant and does not count these changes 
% in the cset matrix. 

% 

if (find(kdofl(chk)=bset))~=Q & kdof2(chk)==’b’ 
base=find(kdoh (chk)=bset); 
kb(base) = kb(base) + lk(chk); 

K_vec(kdofl(chk))=K_vec(kdofl(chk)) + lk(chk); 

else 

if kdofl(chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdofl(chk)]; % matrix and formation of new c-set matrix, 
end 

if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
end 

end 

% Qassical method of reformulatting the [K] matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk)) + Ik(chk); 
if kdof2(chk)~=0 & kdof2(chk)~=’b' 
K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)4cdof2(chk)) - Ik(chk); 
K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)4cdofl(chk)) - Ik(chk); 
K_c(kdof2(chk)Jcdof2(chk))=K_c(kdof2(chk)Jcdof2(chk)) + Ik(chk); 
end 

% Classical method of reformulatting [Kb] matrix. 

% 

ifkdof2(chk)='b’ 

kbb=lk(chk); 

Kb(kdofl(chk),kdofl(chk))=Kb(kdofl(chk),kdofl(chk)) -i- kbb; 
end 

changek=input(’Are there any other changes to the stiffness matrix? l=yes 2=no '); 
end 

changec=inputCWould you like to change your damping matrix? l=yes 2=no '); 
while changec=l 

chc=chc+ 1 ; % counter for determining # columns in mapping matrix. 

flag_c=l; % mapping matrix flag 
lc(chc)=input(’input value of added damper (Ibs-secrin) ’); 

cdofl(chc)=input('input constrained dof where 1st end of added damper is applied '); 
cdof2(chc)=input('input dof where 2nd end of added damper is applied, b for base, or 
0 for ground '); 

% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handl^ when asked where the structure is excited. 

% 

if cdofl(chc)~=bset & cdof2(chc)=='b’ 
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while cdofl(chc)~=bset & cdoG(chc)= 'b' 

dispCConnecting this dof to the base is changing the basic model.’) 
disp('If this is what you desire, go back and include this dof as an excitation 

point.’) 

model_change=input(’Was this change correct? l=yes 2=no ’); 
if model_change=l 

error(’Go back and change your basic model.’) 

else 

cdofl(chc)=input(’input constrained dof where 1st end of added damper is 

applied ’); 

cdof2(chc)=input(’input dof where 2nd end of ^ded damper is applied, b 

for base, or 0 for ground ’); 
end 
end 
end 

% Updates the base damper constant and does not count these changes 
% in the cset matrix. 

% 

if (find(cdofl(chc)~ bset))~=[] & cdof2(chc)=’b’ 
base=find(cdof 1 (chc)==bset); 
cb(base) = cb(base) + lc(chc); 

C_vec(cdofl(chc))=C_vec(cdofl(chc)) + lc(chc); 
else 

if cdofl(chc)~=cset % comparison of change dof with the c-set 

cset=[cset cdofl(chc)]; % matrix and formation of new c-set matrix, 
end 

if cdof2(chc)-'=cset 
cset=[cset cdof2(chc)]; 
end 
end 

% Classical method of reformulatting the [C] matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc)) + lc(chc); 
if cdof2(chc)~=^ & cdof2(chc)~=’b’ 
C_c(cdofl(chc),cdof2(chc))=^_c(cdofl(chc),cdof2(chc)) - lc(chc); 
C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc)) - lc(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 

% Classical method of reformulating [Cb] matrix. 

% 

ifcdof2(chc)=’b’ 

cbb=lc(chc); 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc)) + ebb; 
end 

changec=input(’Are there any other changes to the damping matrix? l=yes 2=no ’); 
end 
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changem=inputCWould you like to change your mass matrix? l=yes 2=no ’); 
while changem— 1 
chm=chm+l; 
flag_m=l; 

lm(chm)=input('input value of added mass (lbs-sec^2/in) '); 
mdof(chm)=dnput(’input constrained dof where mass is added 

% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 

% Classical method of reformulatting the [M] matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + lm(chm); 

changem=input('Are there any other changes to the mass matrix? l=yes 2=no '); 
end 

% Handling of spring added to ground 
% 

ground = find(cset==0); 
cset(ground) = Q; 



% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 

% 

iset=inputCWhat dofs other than those where change occured, are you interested in? ’); 

% Formation of zset (c U b) and eset (i U c U b) matrices. 

% 

zset=[cset bset] ; eset=[iset cset bset] ; 

% Choosing the FRF interested in 
% 

resp_dof=inputCWhat response dof are you interested in? '); 
fff_index=find(eset==resp_dof); 

% Formation of [del_Kc], [del_Cc], [del_Mc], [del_Zc], [del_Zb] 

% & [del_Z] 

% 

del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 

if flag_k=l 
for a=l:chk 
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if (find(kdofl(a)=bset))~=[] & kdof2(a)='b’ 

indexkl=Q; % provides for base changes to not be 

indexk2=D; % included in del_Kc matrix 

else 

indexkl=find(kdofl(a)=cset); % finds where in cset matrix 
indexk2=find(kdof2(a)==cset); % the change in k occurs 
end 

del_Kc(indexkl,indexkl)=del_Kc(indexkl,indexkl) + lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2) - lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 
end 
end 

if flag_c=l 
for a=l:chc 

if (find(cdofl(chc)=bset))~=n & cdof2(a)==’b’ 
indexcl=[]; 
indexc2=Q; 
else 

indexc l=find(cdofl (a)=cset); % finds where in cset matrix 
indexc2=find(cdof2(a)==cset); % the change in C occurs 
end 

del_Cc(indexcl,indexcl)=del_Cc(indexcl,indexcl) + lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2) - lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2,indexcl) - lc(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc24ndexc2) + lc(a); 
end 
end 



if flag_m=l 
for a=l:chm 

indexm=find(mdof(a)==cset); % finds where in cset matrix 

% the change in k occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 
end 
end 

del_Kb=diag(kb); % diag is used in the event >1 base 

excitation 

del_Cb=diag(cb); % is present 

% 

% Time Step: 



start_t = 0.0; 
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%del_t =0.0005; 
%end_t = .04; 

del_t =0.0003; 
end_t = .03; 

%del_t =0.001; 
%end_t = .2; 



% plate times 
% beam times 
% spring times 



time = [start_t:del_t:end_t]; % Time points 
nstep = length(time); % No. Time points 



% 



% Calculate Impulse Response Functions : 
% 



syn_t0 = clock; 

[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi Mmodal zeta 



% 

% Setup and Solve Integral Equation for x2*(t); 

% 

A = zeros(nstep); 

globalA = zeros(length(zset)*size(A,l)); 
count=0; 

col = l+count:nstep4count; 
for acnt = 1 :size(himp, 1 ); 

for icnt_rows = 2 : nstep; 

[wts] = fTrapzWts(icnt_rows); 

A(icnt_rows,l:icnt_rows) = del_t * wts .* fliplr(himp(acnt,l:icnt_rows)); 
end 



row = col(l)4count:col(length(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count = count +nstep; 

if row(length(row)) = size(globalA) 
col = col + nstep; 
count = 0; 
end 



end 



161 



clear A 



% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 
% eset, and these rows and columns are deleted from H_star. This is done 
% because the redundant dof column(s) and row(s) in H_star creates a singular 
% matrix which is not invertable and I&ed cannot be determined. 

% 

for u=l:length(bset) 

locate=find(bset(u)=zset); 
if lengthGocate) > 1 
pull=[puU locate(2)]; 
redundant=[redundant locate(l)]; 
end 
end 

pull_start=pull*nstep-(nstep- 1 ); 

pull_end=pull*nstep; 

delete_dof=[]; 

for v=l:len^(pull) 

delete_dof=[delete_dof pull_start(v):pull_end(v)] ; 
end 

globalA(delete_dof,;)=Q ; 
globalA(: ,delete_dof)=0 ; 



% 

ff = ones(size(globalA,l),l); 

Yo = 1 .0; % Amplitude of base motion 

[Y,Ydot] = fBlastForcing(Yo,time', 'blst', 0); 

tol = le-4; 
dif= 100; 

for icnt = 1:300 



X = -globalA*ff; 



if icnt >= 2 

[ii jj] = max(abs(xlastl - X)); 
dif = max(abs(xlastl(jj) - X(j[j))); 
end 

xlastl = X; 
if dif < tol 

disp('Breaking’); break 
end 

[ff] = Synth_force(X,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 
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end 



syn_tl = clock; 

time_syn_time=etime(syn_t 1 ,syn_tO) 



% 



% Classical Method to Solve for Response: 

% 

clas_tO = clock; 

Xchk = zeros(length(zset),length(time)); 

[wn_c,phi_c,zeta_c,Mmo^_c,phi_norm_c] = modal(M_c,K_c,C_c); 
%zeta_c=0; 



for icnt_modes = 1 : ndof; 

zeta_mode=zeta_c(icnt_modes); 

[mode_irf] = fModalIRF(wn_c(icnt_modes), zeta_mode, time); 
toodal = phi_nonn_c(:,icnt_modes)' * K_vec * Y; 

Xchk = Xchk + phi_norm_c(zset,icnt_modes) * fConvTrap(mode_irf,fmodal',del_t); 
end 

clas_tl = clock; 

clas_syn_time=etime(clas_tl ,clas_tO) 

% Plotting of Transient Response using Classical & Frequency Synthesis Methods 
% 



resp_index=find(zset==resp_dof); 



plot(time,Xchk(resp_index( 1 ),:),’r“',time^((l :nstep)+(resp_index( 1)- 1 )*nstep),'b') 
grid 

title([Transient Time Response at dof ’,int2str(resp_dof)]) 
ylabel('displacement (in)') 
xlabeK'time (sec)') 
legend('Classicar,'Synthesis') 
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tsynode45.m 



% The purpose of this program is to perform the time domain synthesis on a structure and 
% compare the dynamic responses to those calculated using a direct integration method. 

% 

clear 

disp(’Select the file which contains your [M], [K], and [C] matrices' ) 
pause(2) 

[M_K_C,p]=uigetfile(’*.mat',’Load [M], [K], [C]’); 
load (M_K_C) 
dofs=l :ndof; 

C=.00*C; 

[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,C); 

changek=2; changec=2; changem=2; % initializes logic variables 

% to 'no' status 



flag_k=0; flag_c=0; flag_m=0; 



chk=0; chc=0; chm=0; % initializes counter for # of changes to matrices 



cset=0; % initializes cset matrix to zero to avoid null index error in 

% handling the spring to ground if no changes are made to the 
% structure. 

Hey_star=n; HV_c=[]; % (re)initializes these matrices to an empty matrix 

% to avoid matrix size differences each time the 
% model is changed. 

K_c=K; M_c=M; C_c=C; % initializes change matrices which will be used in 

% the classical method, to the original matrices. 

K_vec = zeros(2*size(K_c,l),l); 

C_vec = zeros(2*size(C_c,l),l); 



Kb=zeros(size(K,l)); % initializes the spring, and damper to base matrices 

Cb=zeros(size(C,l)); % to zero. 



b=num2str('b'); 
clear M K C 



% allows b to be entered as a numerical value, but 
% declares b a string variable to used for 
% comparisons in logic statements. 



% Designation where base excitation is located, and the spring constant which 
% connects the substructure to the base which is moving. 
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% 

bset=input('At what dof(s) are your structure excited? '); 
for b_spring=l:length(bset); 

kb(b_spring)=input(fprintf(’input the value of the spring constant which connects dof 
%g to the base ’,bset(b_spring))); 
end 

cb=zeros(l Jength(bset)); % initializes the damper constant which connects the 

% substructure to the base which is moving, to zero. 

% Correcting the classical [K] & pCb] matrices to include the spring element 
% connecting the substructure to the base(s). 

% 

for b_spring=l:length(bset); 

K_c(bset(b_spring),bset(b_spring))=K_c(bset(b_spring),bset(b_spring))+kb(b_spring); 
Kb(bset(b_spring),bset(b_spring))=Kb(bset(b_spring),bset(b_spring)) + kb(b_spring); 
K_vec(size(K_c,l)+bset(b_spring))=K_vec(size(K_c,l)+bset(b_spring)) + kb(b_spring); 
end 

changek=inputCWould you like to change your stiffness matrix? l=yes 2=no ')i 
while changek==l 

chk=chk+l; % counter for determining # columns in mapping matrix. 

flag_k= 1 ; % mapping matrix flag 

lk(chk)=input(’input value of added spring (Ibs^) '); 

kdofl(chk)=input('input constrained dof where 1st end of added spring is applied '); 
kdof2(chk)=input('input dof where 2nd end of added spring is applied, b for base, or 

0 for ground ’); 

% Let's the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handled when asked where the structure is excited. 

% 

if kdofl(chk)~=bset & kdof2(chk)=’b’ 

while kdofl(chk)~=bset & kdof2(chk)=='b' 
dispCConnecting this dof to the base is changing the basic model.’) 
disp(lf this is what you desire, go back and include this dof as an excitation 

point.') 

model_change=inputCWas this change correct? l=yes 2=no ’); 
if model_change=l 

errorCGo back and change your basic model.’) 

else 

kdofl(chk)=input(’input constrained dof where 1st end of added damper is 

applied ’); 

kdof2(chk)=input(’input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground ’); 
end 

end 

end 

% updates the base spring constant and does not count these changes 
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% in the cset noatrix. 

% 

if (find(kdofl(chk)=bset))~=0 & kdof2(chk)=='b’ 
base=find(kdof 1 (chk)=bset); 
kb(base) = kb(base) + lk(chk); 

K_vec(size(K_c,l)+kdofl(chk))=K_vec(size(K_c,l)+kdofl(chk)) + Ik(chk); 
else 

if kdof 1 (chk)~=cset % comparison of change dof with the c-set 

cset=[cset kdofl(chk)]; % matrix and formation of new c-set matrix, 
end 

if kdof2(chk)~=cset 
' cset=[cset kdof2(chk)]; 
end 
end 

% Classical method of reformulatting the [K] matrix. 

% 

K_c(kdofl(chk),kdofl(chk))=K_c(kdofl(chk),kdofl(chk)) + Ik(chk); 
if kdof2(chk)~=0 & kdof2(chk)~='b' 
K_c(kdofl(chk),kdof2(chk))=K_c(kdofl(chk)dcdof2(chk)) - lk(chk); 
K_c(kdof2(chk),kdofl(chk))=K_c(kdof2(chk)Jcdofl(chk)) - lk(chk); 
K_c(kdof2(chk)dcdof2(chk))=K_c(kdof2(chk)dcdof2(chk)) + Ik(chk); 
end 

% Classical method of reformulatting [Kb] matrix. 

% 

ifkdof2(chk)==’b' 

kbb=lk(chk); 

Kb(kdofl(chk),kdofl(chk))=Kb(kdofl(chk),kdofl(chk)) + kbb; 
end 

changek=input(’Are there any other changes to the stiffness matrix? l=yes 2=no *)i 
end 

changec=input(’Would you like to change your damping matrix? l=yes 2=no '); 
while changec=l 

chc=chc-i-l ; % counter for determining # columns in mapping matrix. 

flag_c=l; % mapping matrix flag 
lc(chc)=input('input value of added damper (Ibs-sec/in) ’); 

cdof l(chc)=input(’input constrained dof where 1st end of added damper is applied '); 
cdof2(chc)=input('input dof where 2nd end of added damper is applied, b for base, or 

0 for ground ')J 

% Let’s the user know that he is not allowed to make dof(s) to the base changes which 
% are not included in the bset. A change from dof to base constitutes an addition of 
% excitation points. This should handl^ when asked where the structure is excited- 
% 

if cdofl(chc)~=bset & cdof2(chc)=='b’ 

while cdofl(chc)~=bset & cdof2(chc)=='b’ 



166 



dispCConnecting this dof to the base is changing the basic model.') 
dispC’If this is what you desire, go back and include this dof as an excitation 

point.') 

model_change=inputCWas this change correct? l=yes 2=no '); 
if model_change==l 

error('Go back and change your basic model.') 

else 

cdofl(chc)=input('input constrained dof where 1st end of added damper is 

applied '); 

cdof2(chc)=input('input dof where 2nd end of added damper is applied, b 
for base, or 0 for ground ’); 
end 
end 

end 



% Updates the base damper constant and does not count these changes 
% in the cset matrix. 

% 

if (find(cdofl(chc)=bset))~=[] & cdof2(chc)='b' 
base=find(cdof 1 (chc)=bset) ; 
cb(base) = cb(base) + lc(chc); 

C_vec(size(C_c,l)+cdofl(chc))=C_vec(size(C_c,l)+cdofl(chc)) + lc(chc); 

else 

if cdofl(chc)~=cset % comparison of change dof with the c-set 

cset=[cset cdofl(chc)]; % matrix and formation of new c-set matrix, 
end 

if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 

% Classical method of reformulatting the [C] matrix. 

% 

C_c(cdofl(chc),cdofl(chc))=C_c(cdofl(chc),cdofl(chc)) + lc(chc); 
if cdof2(chc)~=^ & cdof2(chc)~='b' 
C_c(cdofl(chc),cdof2(chc))=Ki:_c(cdofl(chc),cdof2(chc)) - lc(chc); 
C_c(cdof2(chc),cdofl(chc))=C_c(cdof2(chc),cdofl(chc)) - lc(chc); 
C_c(cdof2(chc),cdof2(chc))=C_c(cdof2(chc),cdof2(chc)) + lc(chc); 
end 

% Classical method of reformulating [Cb] matrix. 

% 

if cdof2(chc)=='b' 
cbb=lc(chc); 

Cb(cdofl(chc),cdofl(chc))=Cb(cdofl(chc),cdofl(chc)) + ebb; 
end 

changec=input('Are there any other changes to the damping matrix? l=yes 2=no '); 
end 
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changem=inputCWould you like to change your mass matrix? l=yes 2=no ’); 
while changem=l 
chm=chm+l; 
flag_m=l; 

lm(chm)=input('input value of added mass (lbs-sec'^2;^n) '); 
mdof(chm)=input('input constrained dof where mass is added ’); 

% comparison of change dof with the c-set matrix and formation of 
% new c-set matrix. 

% 

if mdof(chm)~=cset 
cset=[cset mdof(chm)]; 
end 

% Classical method of reformulatting the [M] matrix. 

% 

M_c(mdof(chm),mdof(chm))=M_c(mdof(chm),mdof(chm)) + lm(chm); 

changem=input(’Are there any other changes to the mass matrix? l=yes 2=no '); 
end 

% Handling of spring added to ground 
% 

ground = find(cset=0); 
cset(ground) = Q; 



% Formulation of iset matrix which is the set of dofs other than 
% those in the cset matrix where response information is desired. 

% 

iset=inputCWhat dofs other than those where change occured, are you interested in? ’)5 

% Formation of zset (c U b) and eset (i U c U b) matrices. 

% 

zset=[cset bset] ; eset=[iset cset bset] ; 

% Choosing the FRF interested in 
% 

resp_dof=input('What response dof are you interested in? '); 
frf_index=find(eset==resp_dof); 

% Formation of [del_Kc], [del_Cc], [del_Mc], [del.Zc], [del_Zb] 

% & [del_Z] 

% 

del_Kc=zeros(length(cset)); % intializes the change matrices 
del_Cc=zeros(length(cset)); % to zero and sets their sizes 
del_Mc=zeros(length(cset)); % as [cset x cset] matrix 

if flag_k=l 
for a=l :chk 

if (find(kdofl(a)=bset))~=[] & kdof2(a)=’b’ 
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indexkl=[]; % provides for base changes to not be 

indexk2=Q; % included in del_Kc matrix 

else 

indexkl=fmd(kdofl(a)=cset); % finds where in cset matrix 
indexk2=find(kdof2(a)==cset); % the change in k occurs 
end 

del_Kc(indexkl,indexkl)=del_Kc(indexkl,indexkl) + lk(a); 
del_Kc(indexkl,indexk2)=del_Kc(indexkl,indexk2) - lk(a); 
del_Kc(indexk2,indexkl)=del_Kc(indexk2,indexkl) - lk(a); 
del_Kc(indexk2,indexk2)=del_Kc(indexk2,indexk2) + lk(a); 
end 
end 

if flag_c~l 
for a=l:chc 

if (find(cdofl(chc)=bset))~=D & cdof2(a)=='b' 
indexcl=Q; 
indexc2=G; 

else 

indexcl=find(cdofl(a)==cset); % finds where in cset matrix 
indexc2=find(cdof2(a)=cset); % the change in C occurs 
end 

del_Cc(indexc 1 ,indexc 1 )=del_Cc(indexc 1 .indexc 1 ) + lc(a); 
del_Cc(indexcl,indexc2)=del_Cc(indexcl,indexc2) - lc(a); 
del_Cc(indexc2,indexcl)=del_Cc(indexc2, indexc 1) - lc(a); 
del_Cc(indexc2,indexc2)=del_Cc(indexc24ndexc2) + lc(a); 
end 
end 



if flag_m=l 
for a=l:chm 

indexm=find(mdof(a)==cset); % finds where in cset matrix 

% the change in k occurs 

del_Mc(indexm,indexm)=del_Mc(indexm,indexm) + lm(a); 
end 
end 



del_Kb=diag(kb); % diag is used in the event >1 base excitation 

del_Cb=diag(cb); % is present 



% 

% Time Step: 

start_t = 0.0; 

%del_t = 0.00025; % plate times 

%end_t = .02; 
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del_t =0.0003; 
end_t = .03; 



% beam times 



%del_t =0.001; % spring times 

%end_t = .2; 

time - [start_t:del_t:end_t]; % Time points 
nstep = length(time); % No. Time points 



% 

% Calculate Impulse Response Functions ; 
% 



syn_t0 = clock; 

[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi Mmodal zeta 



% 

% Setup and Solve Integral Equation for x2*(t): 
% 



A = zeros(nstep); 

globalA = zeros(length(zset)*size(A,l)); 
count=0; 

col = l+count:nstep+count; 
for acnt = l:size(himp,l); 

for icnt_rows = 2 : nstep; 

[wts] = fTrapzWts(icnt_rows); 

A(icnt_rows,l;icnt_rows) = del_t * wts .* fliplr(himp(acnt,ldcnt_rows)); 
end 



row = col(l)+count:col(length(col))+count; 
globalA(row,col)=A; 
globalA(col,row)=A; 
count = count +nstep; 

if row(length(row)) == size(globalA) 
col = col + nstep; 
count = 0; 
end 



end 



clear A 

%disp(sprintf('Norm(globalA) = %5.3f ,norm(globalA))) 
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% 

ff = onesOength(zset)*nstep,l); 

Yo = 1.0; % Amplitude of base motion 

[Y,Ydot] = fBlastForcing(Yo,time', ’blst', 0); 

tol = le-4; 
dif=100; 

for icnt = 1:300 

X = -globalA*ff; 



if icnt >= 2 

[ii Jj] = max(abs(xlastl - X)); 
dif = max(abs(xlastl(jj) - X({j))); 
end 

xlastl = X; 
if dif < tol 

disp(’Brealdng’); break 
end 

[ff] = Synth_force(X,Y,Ydot,del_Kc,dd_Cc,del_Mc,kb,cb,nstep,del_t,bset,cset); 



end 

syn_tl = clock; 

time_syn_time=etime(syn_t 1 ,syn_t0) 

% 

% Classical Method to Solve for Response: 

% 

clas_t0 = clock; 

save structure K_c C_c M_c K_vec C_vec Yo 
Xchk0=zeros(2*ndof, 1 ) ; 

[tchk^chk]=ode45(’scnd_to_frst’,start_t,end_t,XchkO); 
clas_tl = clock; 

clas_sy n_time=etime(clas_t 1 ,clas_t0) 



% Plotting of Transient Response using Classical & Frequency Synthesis Methods 
% 

resp_index=find(zset==resp_dof); 
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plot(tchk,Xchk(;,resp_dof(l)),’r— ',time,X((l:nstep)+(resp_index(l)-l)*nstep),'b') 
grid 

titie([Transient Time Response at dof ',int2str(resp_dof)l) 
ylabel('displacement (in)') 
xlabeK'time (sec)') 
legend('ClassicaJ', 'Synthesis') 
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tsynmain.m 



% This program is designed to calculate the new transient responses of a system after 
% changes in the mass, stiffness, and damping matrices have been made. 

% The new responses will be calculated using the synthesis method in the time domain. 
% 

clear 

load change 

% Designation of the frequency range and step size 
% 

% initial step size end 

tjnput = [ 0.0 0.0004 .07 ]; 

time = t_input(l):t_input(2):t_input(3); 

nstep = length(time); % No. Time points 

[wn,phi,zeta,Mmodal,phi_norm] = modal(M,K,C); 
clearMKC 

syn_t0 = clock; 

[himp] = impmodal(wn,phi,zeta,Mmodal,time,zset); 
clear wn phi zeta Mmodal 

[globalA] = buildA(himp,bset,zset,nstep,t_input(2)); 

Yo = 1.0; % Amplitude of forcing function 

[Y,Ydot] = fBlastForcing(Yo,time', 'blst', 0); 

[X_star,icnt] = timesynth(globalA,zset, Y, Y dot,del_Kc,del_Cc,del_Mcdcb,cb ,. . . 

nstep,t_input(2),bset,cset); 



syn_tl = clock; 

time_syn_time=etime(syn_t 1 ,syn_t0) 

% Plotting of Transient Response using Classical & Frequency Synthesis Methods 
% 

resp_dof=input('What response dof are you interested in? ')» 

[X_star_desired] = time_sift(zset 4 'esp_dof,X_star,nstep); 

plot(time,X_star_desired,'b') 

grid 

title(['Synthesized Transient Time Response at dof ’,int2str(resp_dof)]) 
ylabel(’displacement (in)') 
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xlabel('time (sec)') 



impmodal.m 



function [himp] = impmodal(wn,phi,zeta,Mmodal,t,imp_dof) 

% This function is designed to calculate the IRF using mode shapes and 
% summing the IRF over a number of modes. 

% 

ncol=(length(imp_dof)*(length(imp_dof)+l))/2; 

% Calculation of impulse response function 
% 

himp=zeros(ncol,length(t)); % Initialization to zero, and sizing 

% of the freq x lower triangle matrix. 

nmodes=size(phi); b=0; % Defining the number of modes that exist and 

% initializing b which keeps track of the 
% of modes 

for b=l:nmodes 

% Elimination of unnessecary elements and rearrangement of original 
% mode shapes {phi) to form (phi_red) and [num_ied] which is now an 
% imp_dof X imp_dof matrix. 

% 

num_red=phi(imp_dof,b)*phi(imp_dof,b)'; 

if wn(b) > le-3 % Then elastic mode 

wd = wn(b) * sqrt(l - zeta(b)'^2); 

mode_irf = exp(-zeta(b)*wn(b)*t) .* sin(wd * t) / (Mmodal(b)*wd); 

elseif wn(b) <= le-3 % Rigid body mode 
mode_irf = t/Mmodal(b); 
end 

% Changes the [num_red] matrix from an imp_dof x imp_dof matrix to 
% a 1 X lower triangle vector. 

% 

H_lowtri=tril(num_red); % pulls out the lower triangle 

% and places zeros everywhere else 

symtry = find(H_lowtri==0); % finds zero values in the 

% lower triangle matrix 

H_lowtri(symtry) = []; % deletes all values of zero and 
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% turns the remaining elements 
% into a 1 X lengthOower triangle) 
% vector 



for cnt=l :length(H_lowtri) 

himp_m<xie(cnt,:)=H_lowtri(cnt)*mode_irf; 

end 

himp = himp + himp_mode; 
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buildA.m 



function [globalA] = buildA(himp,bset,zset,nstep,del_t) 

% The purpose of this program is to build the trapezoidal quadrature matrices to be used 
% for the solution of the VIDE 
% 



A = zeros (nstep); 

globalA = zeros(length(zset)*size(A,D); 
count=0; 

col = l+count:nstep+count; 
for acnt = l:size(h^p,l); 

for icnt_rows = 2 : nstep; 

[wts] = fTrapzWts(icnt_rows); 

A(icnt_rows,l:icnt_rows) = del_t * wts .* fliplr(himp(acnt,i:icnt_rows)); 
end 

row = col(l)+count:col(length(col))+count; 
globa^(row,col)=A; 
globaiA(col,row)=A; 
count = count +nstep; 

if row(length(row)) = size(globalA) 
col = col + nstep; 
count = 0; 
end 

end 

% Checks to see if there are any redundant dofs in the eset due to changes and 
% bset occurring at the same dofs. If there is redundancy, it is located in 
% eset, and these rows and columns are deleted from H_star. 

% 

for u=l:length(bset) 

locate=find(bset(u)=zset); 
if length(locate) > 1 
pull=[pull locate(2)]; 
redundant=[redund^t locate(l)]; 
end 
end 

pull_start=pull*nstep-(nstep- 1 ); 

pull_end=puU*nstep; 

delete_dof=D; 

for v=l :length(pull) 

delete_dof=[delete_dof pull_start(v) :pull_end(v)] ; 
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end 

globalA(delete_dof,:)=[]; 

globalA(:,delete_dof)=D; 
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fBlastForcing.m 



function [f_of_t,fdot] = fBlastForcing(Fo,time, type, plotit); 

% 

% Usage: [f_of_t,fdot] = fBlastForcing(Fo,time, type, plotit); 
% 

% Choices: sine blst step 
% 

% type = 'step' STRING Variable 

% 

% 

% If use ’sine’, fdot also returned. 

% 

% This function returns a forcing function which is 
% a "blast" function. 

% 

% F(t) = Fo * ( exp(-at) - exp(-bt) ) 

% 

% where a and b are constants which shape the blast, 

% and Fo is the amplitude of the blast. 

% 

% The variable "plotit" is a switch which if set = 1 will 
% cause the f(t) to be plotted, if set to anything else 
% will not plot 
% 

% 

% Choices: sine blst step 
% type = 'step'; 

if type = 'blst'; 

% dispC Blast forcing used...') 
a = 100.0; 
b = 300.0; 

f_of_t = Fo * ( exp(-a*time) - exp(-b*time) ); 
fdot = Fo * ( -a*exp(-a*time) + b*exp(-b*time) ); 



elseif type == 'step'; 



% 



dispC Step forcing used...') 
f_of_t = Fo * ones(size(time)); 
fdot = []; 



elseif type == 'sine'; 
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% dispC Sine forcing used...’) 

W = 5; % Hertz 

= Fo * sin(2*pi*W*time); 
fdot = Fo * (2*pi*^*cos(2*pi*W*time); 



end; 



if plotit = 1; 
figure(gcf+l) 

if type = 'sine'; 

plot(tirne,f_of_t,tinie,fdot);grid 

else 

plot(time,f_of_t);grid 
end 
pause 
% clf 
end 

% End function. 
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timesynth.m 



function [X_star,icnt] = timesynth(globalA,zset,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,... 

nstep,del_t,bset,cset) 

% This function is designed to calculate the new dynamic response using the synthesis 
% method. 

% 

ff = ones(size(globalA,l),l); 

tol = le-2; 

dif=100; 

for icnt = 1:300 

X_star = -globalA*ff; 

if icnt >= 2 

[ii,jj] = max(abs(xlastl - X_star)); 
dif = max(abs(xlastl(jj) - X_star(ij))); 
end 

if icnt>=300 
ii 



end 

xlastl = X_star; 
if dif < tol 

% disp(’Breaking’); 
break 
end 

% icnt 
[ff] = 

Synth_force(X_star,Y,Ydot,del_Kc,del_Cc,dd_Mc,kb,cb,nstep,del_t,bset,cset); 

end 
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time sift.m 



function [X_star_desired] = time_sift(zset4‘esp_dof^_star,nstep) 

% The purpose of this program is to sort through the synthesized transient 
% responses and return the one in which the user is interested in. 

% 

resp_index=find(zset==resp_dof); 

X_star_desired = X_star((l:nstep)+(resp_index(l)-l)*nstep); 
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APPENDIX D. OPTIMIZATION COMPUTER CODES 



The following is a list and a brief description of the main MATLAB 

computer codes that were written in order to perform the optimization of a 

shock and vibration isolation system. 

• isomain.m - the main program which loads the baseline structure to be 
optimized, and allows for the designation of the design variables. Also, if 
necessary, calls various modularized functions such as modal.m, 
frfmodal.m, impmodal.m, fBlastForcing.m, buildA.m. It then allows for 
the optimization inputs, the use of the MATLAB Optimization function 
constr.m, and the post processing. 

• modal.m - solves for the structure's natural frequencies, mode shapes and 
other modal matrices and vectors. 

• frfmodal.m - calculates the structure's FRFs. 

• impmodal.m - calculates the structure's IRFs. 

• buildA.m - builds the quadrature matrices. 

• fBlastForcing.m - inputs the base excitation as a fimction of time. 

• isosub.m - the subroutine program which is called by the MATLAB 
Optimization function constr.m. This is where the optimization iterations 
occur, the design variables change, and the system's responses are 
calculated. This program calls the fimctions delta.m, frfsynth.m. 
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statdelta.m, statsynth.m and timesynth.m. It then calculates the 
normalized constraint values. 

• delta, m - forms the modification matrices which will be used to form 
impedance matrices or determine coupling forces. 

• frfsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized FRFs. 

• statdelta.m - forms the modification matrices which will be used to form 
impedance matrices. 

• statsynth.m - performs the synthesis on the baseline structure and returns 
the synthesized static displacements. 

• timesynth.m - performs the synthesis on the baseline structure and 
returns the synthesized transient responses. 

The full codes are contained on subsequent pages. Any codes that were 

previously presented will not be repeated. 
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isomain.m 



clear 

clear global 
diary isomain.dia 
tic 



global kdofl kdof2 cdofl cdof2 mdof iset cset bset zset eset kbO cbO hee omega ... 
opt_omega excite resp_dof inp_dof resp_index stat_omega stat_hee iter ... 
del_vars del_f H_desired_0 del_H_desired globalA Y Ydot t_input nstep 



% 

% Loading of Structure to be Optimized 

% 

% 



dispCAre there already reduced FRF & IRF matrices in the correct format available,') 
dispC’or does the presynthesized FRF & IRF matrices need to be generated using the dof ’) 
dispC'locations where you desire to change the structure? ') 

FRF_IRF=input(’l=reduced FRF/IRF already exists 2=generate reduced FRF/IRF '); 



% Protects against user making an eiror in choosing how FRF/IRF is obtained. 

% 

if FRF_IRF~=[1 2] 
while FRF_IRF~=[1 2] 

dispCEixor in choosing how FRF/IRF is obtained. Choose 1 or 2.’) 
FRF_IRF=input(’l=reduced FRF/IRF already exists 2=generate reduced FRF/ERF '); 
end 
end 



ifFRF_IRF=l 

dispCSelect the file which contains your presynthesized FRF/IRFs, a corresponding 
frequency ’) 

dispCand time vectors, and the vector of all the dofs that will be in your eset vector.’) 
pause(2) 

[hee_dofs_omega,p]=uigetfile(’*.mat','Load FRF/IRF’) 
load (hee_dofs_omega) 

disp(’Also select the file which contains your presynthesized FRFs with zero damping ’) 

dispCfor the purpose of calculating static displacement.’) 

pause(2) 

[stat_hee,p]=uigetfile(’*.mat’,’Load Static FRF’) 
load (stat_hee) 

else 

dispC’Select the file which contains your [M], [K], and [C] matrices’ ) 
pause(2) 
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[M_K_C,p]=uigetfile(’*.mat’,’Load [M], [K], [C]'); 
load (M_K_C) 
dofs=l:ndof; 
kcO=K(109,109); 
end 



% 

% Designation of Optimization and Synthesis DOFs 

% 

% 

kdofl=D; kdof2=n; 

cdofl=Q; cdof2=Q; 

mdof=D; 



% initalizes the matrices which keep 
% track of the dofs where changes occur 



b=num2str('b'); % allows b to be entered as a numerical value for the 

dof 



used for 



% choices; but declares b a string variable to 
% comparisons in logic statements. 



kdofs=input(’Input dof pairs for stiffness changes for optimization '); 
if kdofs— [] 
kdofl=[]; kdof2=D; 
else 

kdof l=kdofs(:, 1 ); kdof2=kdofs(;,2); 

end 



cdofs=inputCInput dof pairs for damping changes for optimization ’); 
if cdofs=[] 
cdofl=[]; cdof2=Q; 
else 

cdofl=cdofs(:,l); cdof2=cdofs(:^); 

end 



mdof=input(’Input dofs for mass changes for optimization ’); 



^ afe ^ 3{e ^ 3{e )(e 3{e 3)e a{e ^ 3fe a|e 3|e a|e 3{e 3)c 3(e ^ :)c :|e :jc 3|e ^ ^ :|c afe :)c 3|e 3|e a|e :}e 3|e 9|e 9(e 9(c ^ 

% 

% Formation of Partitioning and Arrangement Vectors 

% 

% 



excite=inputCHow is the structure excited? l=system 2=base '); 

% Formation of (bset} 

% 

cset=0; bset=[]; % initializes cset matrix to zero to avoid 

% null index error if no changes are made to 
% the stmeture and bset to empty matrix to 
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% prevent error in the event this is not a 
% base excitation and bset does not exist 

kb0=0; cb0=0; % initializes kb and cb to zero to prevent error in 

% the event this is not a base excitation and kb & cb does 
% not exist. 



if excite=2 

% Designation where base excitation is located. 

% 

bset=input(’At what dof(s) are your structure excited? '); 



% Designation of the spring and damper constants which connects the substructure to the 
% base which is moving.and protection against user making an error in choosing the 
% value of the constant(s). 



% 

for b_spring=l:length(bset) 

kbO(b_spring)=input(fprintf('input the value of the spring constant which connects 

dof %g to the base ’,bset(b_spring))); 

kb0_zero=find(kb0==0); 
if length(kbO)~=b_spring I kbO_zero~=Q 

while length(kbO)~=b_spring I kbO_zero~=[] 

fprintf('You made an error in entering the value for dof %g base excitation 
spring constant.\n',bset(b_spring)) 
kb0(b_spring)=0; 

kbO^_spring)=input(^rintf(’Re-input the value of the spring constant which 
connects dof %g to the base ',bset(b_spring))); 
kb0_zero=find(kb0=0); 
end 
end 
end 

cbO=zeros(l ,length(bset)); % initializes the damper constant which connects the 

% substructure to the base which is moving, to zero, 
end 



% Formation of {cset} vector 
% 

if kdofs~=n 
for chk= 1 ;length(kdof 1 ) 

% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
% Does not count changes to the base spring constant in the cset matrix. 

% 

if excite==2 & (find(kdofl(chk)=bset))~=[] & kdof2(chk)='b' 
cset=cset; 

else 

if kdof l(chk)~=cset 
cset=[cset kdofl(chk)]; 
end 

if kdof2(chk)~=cset 
cset=[cset kdof2(chk)]; 
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end 



end 

end 

end 

if cdofs~=[] 

for chc=l :length(cdofl) 

% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
% Does not count changes to the base damper constant in the cset matrix. 

% 

if excite=2 & (find(cdofl(chc)=bset))~=[] & cdof2(chc)=='b' 
cset=cset; 

else 

if cdofl(chc)~=cset 
cset=[cset cdofl(chc)]; 
end 

if cdof2(chc)~=cset 
cset=[cset cdof2(chc)]; 
end 
end 
end 
end 



ifmdof~=[] 

for chm=l:length(mdof) 

% Comparison of change dof with the c-set matrix and formation of new c-set matrix. 
% 

if mdof(chm)~=cset 

cset=[cset mdof(chm)]; % and formation of new c-set matrix, 

end 
end 
end 

ground = find(cset==0); % Handling of spring added to ground to 

eliminate 

cset(ground) = []; % zero from cset 

% Formation of {iset} vector 
% 

ifFRF_IRF=l 

dispCSince you inputed your own FRF/IRF matrix, your iset is chosen as all dofs 
remaining ') 

dispCin your eset/dof vector after extracting the cset vector.’) 
iset=dofs; 
for a=l:length(cset) 
extract(a)=fmd(cset(a)=dofs); 
end 

iset(extract)=n; 

else 
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% iset=inputCWhat dofs other than those where change occured, are you interested in? 

iset=[]; 

end 



% Formation of zset (c U b) and eset (i U c U b) matrices. 
% 

zset=[cset bset]; eset=[iset cset bset]; 



^ :jc :jc :fc 4c % ?ic :jc ?ic ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ :jc ;jc :jc ?ic :jc :jc :jc % % 4c % 4c :jc 4c ?ic:fc :fc :fc :|c :jc :ic :jc :ic :fc :jc 

% 

% Formulation of original [H], [Himp], for Freq & Time Synthesis and Static 
% Displacement 
% 



ifFRF_IRF=2 

% Designation of the frequency and time range and step size 
% 

% initial step size end 

freq = [ .01 10 le3+.01 ]; 

t_input = [ 0.0 0.001 .05 ]; 

omega=freq(l):freq(2):freq(3); 

time = t_input(l):t_input(2):t_input(3); 

nstep = length(time); % No. Time points 

[wn,phi,zeta,Mmodal,phi_norm,Cmodal] = modal(M,K,C); 
[hee] = frfmodal(wn,phi,zeta,Mmodal,omega,eset); 

[himp] = impmodal(wn,phi,zeta,Mmo^,time,zset); 

% Information for static displacement synthesis 
% 

stat_zeta = 0*zeta; 
stat_omega = . 1 :. 1 :.4; 

[stat_hee] = frfmodal(wn,phi,stat_zeta,Mmodal,stat_omega,eset); 
end 

clear M K C 



^ 4c4e 4e 4 c 4 c 4c 4e 4c 4e 4c 4c 4c 4e 4: ** 4c 4c 4c 4c 4c ************ 4c 4e 4:** * 4c 4c 4c 4c 4c 4c ** 4c 4c 4i ** 4c 4c 4c 4c 4c 4c 4c 4c * 4c 4c 4c 4c 4c 4c 4c 4c 4c 



% 

% Time Synthesis 
% 



% Builds the matrix A which uses trapezoidal weights to integrate the IRFs 
% 

[globalA] = buildA(himp,bset,zset,nstep,t_input(2)); 
clear himp 
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% Fomas the base excitation vectors of displacement and velocity 

% 

Yo = 1.0; % Amplitude of forcing function 

[Y,Ydot] = fBlastForcing(Yo,time’, 'blst', 0); 



% 

% The Optimization 
% 



% Designates the frequency over which the optimization will be conducted 
% 

if FRF_IRP=2 

dispCSince you generated the FRF/ERFs through this optimization program, your FRF’) 
dispCfrequencies and IRF times are automaticdly chosen as the freq. and time range ’) 
dispCfor your optimization.') 
opt_omega=omega; 
else 

opt_freql=input('What is the initial frequency over which to evaluate this problem’); 
opt_freq2=input(’What is the final frequency over which to evaluate this problem’); 
dispCThe frequency increment is automatic^ly chosen to match the FRF freq. 
increment’) 

opt_omega=opt_freq 1 :omega(2)-omega( 1 ):opt_freq2; 
end 



% Selecting the FRF that is to be minimized. 

% 

resp_dof=input(’What response dof are you interested in? ’); 

if excite=2 
inp_dof=D; 
else 

inp_dof=input('What input dof are you interested in? ’); 
end 

resp_index=find(zset==resp_dof); 

% Setting the number of of optimization variables, their initial values, and 
% their upper and lower bounds. 

% 

numoptvar=input(’How many optimization variables do you have? ’); 
del0=[l 1 1 1]; 

% del_kp del_kc del_mp 

% del_kb deljcc del_cb del_cc 

dellb= [ -4 -4 0 0 ]; 

delub= [ 5 5 5 5 ]; 
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% Initializing the matirices which will keep track of how the design variables 

% and objective function changes for every optimization iteration for the 

% purpose of plotting. 

% 

del_vars = []; del_f = []; del_H_desired = □; 

iter=0; % Sets the number of iterations counter at zero 

% Calling of constr.m routine to perform optimization 

% on subroutine iso_sub.m 

% 



options(l)=l; 

options(14)=200; 

toe 



tic 

[del]=constr('isosub',delO, options, ddlb.delub); 



^ ;{e ;fc ;ic sic ;{c 4c ;{e % :{c ;{c ;|c ;fc 4c 3k 3k 3k ^ 3|c % % :|c 3{c ;|c 3k 3fc 4c 3k 3k :fe He 4c 3fc 3fc 9k ^ ^ 3k % 



% OUTPUT 

% 

load opt_data 

final_del_kb = del_vars(length(del_vars),l); 
final_del_kc = del_vars(length(del_vars),2); 
final_del_cb = del_vars(length(del_vars),3); 
final_del_cc = del_vars0ength(del_vars),4); 

disp(The recommended change in base isolator stiffness GbAn) is — > ’),final_del_kb 

disp(The recommended change in computer isolator stiffness GbAn) is — > ’),fmal_del_kc 

disp(The recommended change in base isolator damping (IbAn/s) is — > '),final_del_cb 

dispCThe recommended change in computer isolator damping Gb/in/s) is — > 

’),final_del_cc 

dispCConstraints:’) 

disp(g) 

iterations=l titer, 

% Plot of FRF to observe how it changes during the optimization 
% 

figure(l) 

semilogy(omega,abs(del_H_desired(l, r',omega,abs(del_H_desired(2,:)),'m') 

title('Optimized Frequency Response Function [H*]’) 

ylabelCFRF Amplitude (unitless)') 

xlabelf'Frequency (rad/s)’) 

legend('Before’,'After’) 

figure(2) 
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subplot(2,l,l) 

plot(iteradons,del_vars(:,l),'-r', iterations, del_vars(:, 2), 'm') 
title('Change in Variables vs. Iterations’) 
ylabel (Isolator Stiffness (Ib/in)') 
legend(’deLkb','del_kc’) 
subplot(2,l,2) 

plot(iterations,del_vars(:,3),’— r’,iterations,del_vars(:,4),'m’) 

titleC’Change in Variables vs. Iterations’) 

xlabel(’# of iterations’) 

ylabel(lsolator Damping (Ib-s^n)’) 

legend(’del_cb’,’del_cc’) 

figure(3) 

subplot(2,l,l) 

plot(iterations,del_vars(:, l)+kbO(l),’-r’, iterations, del_vars(;,2)+kc0,’m’) 
title(’Change in Isolation Constants vs. Iterations’) 
ylabel(’Isolator Stiffness (Ib/in)’) 
legend(’k_base’,’k_computer’) 
subplot(2,l,2) 

plot(iterations,del_vars(: ,3),’— r’,iterations,del_vars( : ,4),’m’) 
title(’Change in Isolation Constants vs. Iterations’) 
xlabel(’# of iterations’) 
ylabel(’Isolator Damping (Ib-s/in)’) 
legend(’c_base’,’c_computer’) 

figure(4) 

semilogy(iterations,del_f,’r’) 

title(’Change in Objective Function vs. Iterations’) 

ylabel(’Magnitude’) 

xlabel(’# of iterations’) 



figure(5) 

plot(time,X_star((l:nstep)+(resp_index(l)-l)*nstep),’b’) 

grid 

title([’Synthesized Transient Time Response at dof ’,int2str(resp_dof)]) 
ylabel(’displacement (in)’) 
xlabel(’time (sec)’) 

figure(6) 

plot(time,Xaccel_star/386.4,’m’) 

grid 

title([’Synthesized Transient Time Response at dof ’,int2str(resp_dof)]) 
ylabel(’acceleration (g"s)’) 
xlabel(’time (sec)’) 



^ 3k 3k He ^ sic ^ :f: :je :jc :f: :f; :f: 3): :f: 3): 3k 4: :{c :je :f: :fe 3^ 



toe 

diary off 

% end Iso_main.m 
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isosub. m 



function [f,g]=iso_sub(del) 

global kdofl kdofZ cdofl cdof2 mdof iset cset bset zset eset kbO cbO hee omega... 

opt_omega excite resp_dof inp_dof resp_index stat_omega stat_hee iter del_vars... 
del_f H_desired_0 del_H_desired globalA Y Ydot t_input nstep 

iter=itCT +1; % counts number of iterations 



% Initialization and assigning of the optimization variables to changes which are made to 
% the structure. Spring and damper constants are scaled to even out the domain. 

% 

lk=0; lc=0; lm=0; 

del_kb=le3*del(l); del_kc=le3*del(2); 

del_cb=l*del(3); del_cc=l*del(4); 



lk(l:4)=del_kb*ones(l,4); lk(5)=del_kc; 

lk(6:9)=- 1 e3*ones(l ,4); 

lc( 1 :4)=del_cb*ones( 1 ,4) ; lc(5)=del_cc; 



^ :{c :ic :{c 9k sk 3k ^ ^ % 4^ 9k Ne sk sk ^ sic :ic :}c 3k ^ 3k 9k Ne sic 3k 3k 9k 3k ^ 9k ^ 3k ^ 9k 3k 9k 9k ^ 9k ^ 3k 3|c :k 3k 3k 9ic ^ 3k 3k 9k 3^ 

% 

% Formulation of FRF and Calculation of mean square acceleration 
% 



[del_Kc,del_Cc,del_Mcdcb,cb] = 

delta(cset, bset, kdofl ,kdof2,lk,cdof 1 ,cdof2,lc,mdof,lm,kb0,cb0); 

[h_star] = frfsynth(hee,kb,cb,omega,excite,eset,iset,bset,del_Kc,del_Cc,del_Mc); 
[H_desired] = fif_sift(eset,resp_dof,inp_dof,excite,h_star); 

cut_fireq=[find(omega<opt_omega( 1 )) find(omega>opt_omega(length(opt_omega)))] ; 
H_desired(cut_freq)=[] ; 

H_objective=(abs(H_desired)).'^2; 

[f]=simpl_3(H_objective,opt_omega(2)-opt_omega(l)) 



5k 3k 9k 3k 3k 3k 3k 9k 9k 3k 9k 3k 3k 3k 3k 9k 9k 3k 9k 3k 3k 9k 3k 9k 9k 9k 9k 9k 9k Jk Jk 3k 3k 9k 9k 3k 3k 3jc 3k 9k 3k 3k sk 9k sk 3k 9k 3k 3k 3k 9k 3k 3k 3k 9k 9k 3k 9k 9k 3k 9k 9k 9k 9k 9k 9k 9k 3k 9k 9k 



% 

% Plot Generation 
% 



% Saving a matrix of the optimization variables, objective 
% plot how they changed during the optimizatiion. 



% 



function, and FRF in order to 
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if iter=l 

del_vars = [del_kb del_kc del_cb del_cc]; 
H_desired_0 = H_desired; 
deLf=[f]; 
else 

del_vars = [del_vars;del_kb del_kc del_cb del_cc]; 
del_H_desired = [H_desired_0; H_desired]; 

del_f = [del_f;f]; 
end 



clear hee H_desired_0 H_desired 



^4c*:jc4c4c*4c:ic4c4c**4c:ic**4c4c4c4c:ic*4c**:jc:|c:ic:ic4c4c4c:jc4c:ic4c4c:ic4c*4c4c*4c:ic4c*%*%4c******3|c****3ic* ****** 
% 

% Calculation of plate and computer static defection 
% 



[stat_del_Kc,stat_del_Mc,stat_kb] = statdelta(cset,bset,kdof 1 Jcdof2,lk,mdof,lmJcbO); 

[z_stat,Fred,Kred,Mred,HstarO_dp] = statsynth(stat_hee,stat_kb,stat_omega 

eset,iset,bset, stat_del_Kc,stat_del_Mc) ; 



% 



% Calculation of dynamic response due to shock 



[X_star,icnt] = timesynth(globalA,zset,Y,Ydot,del_Kc,del_Cc,del_Mc,kb,cb,... 

nstep,t_input(2),bset,cset); 



fprintfCicnt = %g',icnt) 
%dispCFinish time synthesis') 



% 



% Calculation of normalized constraint values 



maxzp_stat=-.08; maxzc_stat=-.03; 

maxkp_stretch=l; maxkc_stretch=l; 

maxzc_accel=30* 3 86.4; 



zp_stat=z_stat([l 3:6]); % puUs out the static displacements for the plate 

zp_stat=-(max(abs(zp_stat))); % and determines where the maximum displacement 

% is min command used because static displacement 
% is negative 



zc_stat=-(abs(z_stat(2)) - abs(z_stat(l))); 
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% calculation of the relative static 
% displacement between the plate and 
% computer 



% calculation of the maximum relative dynamic displacements between the plate and the 
% base and the computer and the plate 
% 

kp_stretch = X_star(2*nstep+l:6*nstep) - [Y;Y;Y;Y]; 
kp_stretch = max(abs(kp_stretch)); 
kc_stretch = X_star(nstep+l:2*nstep) - X_star(l:nstep); 
kc_stretch = max(abs(kc_stretch)); 

% calculation of the maximum absolute acceleration of the computer 
% 

Xaccel_star = fddot(X_star((l:nstep)+(resp_index(l)-l)*nstep),t_input(2)); 
zc_accel = max(abs(Xaccel_star)); 

g(l)=zp_stat/maxzp_stat - 1; 
g(2)=zc_stat/maxzc_stat - 1; 
g(3)=kp_stretch/maxkp_stretch - 1; 
g(4)=kc_stretch/maxkc_stretch - 1; 
g(5)=2c_accel/maxzc_accel- 1 ; 

save opt_data del_vars del_f iter del_H_desired g z_stat kp_stretch kc_stretch zc_accel ... 
X_star XacceLstar 
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